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Abstract 



In this thesis I study the dynamics of a collapsing scalar field coupled to Einstein's 
equations. In this model, evolution of initial data leads to one of two possible endstates: 
formation of a black hole or dispersion to flat space. At the threshold between black 
hole formation and dispersion the behavior of the system is characterized by so-called 
critical phenomena: scaling, self-similarity and universality. These features of critical 
collapse are numerically investigated from both local and global points of view. 

On the one hand, only a small region of spacetime close to the origin is relevant 
for the dynamics of critical collapse. On the other hand, it is also possible for a distant 
observer to analyze the radiation signal emitted by the collapsing matter field. In the 
framework of characteristic evolution, such observers can be modelled by employing 
radial compactification on outgoing null cones, so that the numerical grid extends to 
future null infinity. One may then extract global properties such as the Bondi mass and 
the news function. 

We study the threshold behavior by numerical evolution of one parameter families of 
initial data. The parameter is fine-tuned to the threshold via bisection. In the evolution 
of such near-critical data, the solution approaches the self-similar critical solution for 
some time. We find that the critical exponent that characterizes the discretely self-similar 
solution can be extracted both locally, in the self-similar region and globally, e.g. from 
the news function. In this sense, self-similarity is observable from future null infinity. 

For late times in near-critical evolutions, we see a residual mass concentrated outside 
of the self-similar region and conjecture that it originates from backscattering of outgoing 
radiation during the collapse. The fate of this mass is unclear. If, in supercritical 
evolutions, this mass were to fall into the black hole, the black hole mass would be 
finite, no matter how fine-timed the initial data. 

For subcritical evolutions we have numerically analyzed the exponents of power law 
tails and have found agreement with analytical calculations for radiation along null 
infinity and along timelike lines. We argue that for astrophysical observers the relevant 
falloff rate is that of future null infinity. 

We have also investigated the behavior of quasinormal modes (QNM) in near-critical 
evolutions. Although the perturbation theory underlying QNMs requires a fixed black 



iii 



hole background, we have found a surprising correlation between the radiation signal 
with the period of the first QNM determined by the time-dependent Bondi mass. In 
this context, we have also been able to verify a stationarity criterion based on QNMs for 
the Vaidya metric, which models the time-dependent Schwarzschild backgroimd. 
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Chapter 



1 



Introduction 



In the context of general relativity, consider the collapse of a spherical shell of matter 
under its own weight. The dynamics of this process, as modelled by the coupled Einstein 
and matter field equations, can be understood intuitively in terms of the competition 
between gravitational attraction and repulsive internal forces (due, for instance, to 
kinetic energy or pressure). Typically, such an isolated system ends up in one of 
three distinct states. If the initial configuration is dilute, then the repulsive forces 
will dominate and the collapsing matter will implode through the center and disperse, 
leaving flat space behind. If, on the other hand, the density of the initial configuration is 
sufficiently large, some fraction of the initial mass will form a black hole. In some matter 
models it is also possible to form stable stars, but, for the sake of simplicity, we will 
disregard this possibility in the following. Critical gravitational collapse occurs when 
the attracting and repulsive forces governing the dynamics of the collapse process are 
almost in balance, or, in other words, the initial configuration is near the threshold of 
black hole formation. 

Critical phenomena in gravitational collapse have been originally discovered in the 
seminal numerical investigations of scalar field collapse by Choptuik [dT92a, Cho93]. 
Using sophisticated numerical techniques, Choptuik investigated the threshold of black 
hole formation for the self-gravitating massless scalar field in spherical symmetry. He 
evolved one-parameter families of initial data that interpolate between black hole for- 
mation and dispersion and fine-tuned the initial data parameter, p, through a bisection 
search to its critical value, p*, where a black hole is just formed. Choptuik was able 
to give convincing evidence that black holes of arbitrarily small mass can be created. 
Moreover, he discovered the following surprising phenomena: The black hole mass 
depends on the initial data parameter via a simple power law 



for p > p*. All near-critical solutions approach a discretely self-similar solution at 
intermediate times. This so-called "critical solution" or "Choptuon" is characterized by 
a constant A. These phenomena and the "critical exponents" y and A are independent 
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of the family of initial data. Therefore, critical phenomena are universal within a given 
model. In the dynamical systems picture, the phase space (or space of initial data) of 
this system is divided into basins of attraction, with black holes and Minkowski space 
as attractors. 

Hamade and Stewart [HS96] have found numerical evidence, that the critical solu- 
tion contains a naked singularity which can be seen at future null infinity. The problem 
has also been studied extensively from an analytic point of view by Christodoulou 
[Chr86, Chr87, Chr91, Chr94]. In particular, he was able to prove that the space 
of regular initial data that lead to naked singularities has measure zero [ ], so 
that the appearance of the naked singularity is non-generic. Similar critical solutions 

- exhibiting (continuous or discrete) self-similarity - have also been found for sev- 
eral other types of matter fields, and have been constructed directly in several cases 
[Gun95, Gun97b, Gun97a, LTHA02, LecOl]. 

In this thesis we present further numerical studies of spherically symmetric scalar 
field critical collapse. We extend previous investigations by focussing on global aspects 
of this problem, and, for the most part, use a compactified evolution scheme which 
includes null infinity on our numerical grid. The motivation is twofold: The main goal 
of our investigation was to gain an understanding of local versus global issues in critical 
collapse. In particular, we try to address questions like: What is the role of asymptotic 
flatness for critical collapse (e.g. the critical solution, the "Choptuon" is self-similar, 
and thus not asymptotically flat)? How would hypothetical detectors of radiation 
observe the dynamics close to criticality? How can we understand the way null infinity 
approximates observers at large distances in this simple but nontrivial setup? The 
second motivation is to test numerical algorithms which are based on compactification 
methods in a situation that is very demanding on accuracy. We will argue that at least 
in the model considered here, global methods do not cause a significant penalty in 
accuracy, but simplify the interpretation of certain results. 

In the current work, we refer to critical collapse phenomena as "critical collapse 
at the threshold of apparent horizon formation" to avoid possible misunderstandings, 
since critical collapse is essentially a quasilocal phenomenon and the standard definition 
of black holes is based on global concepts (see textbooks like e.g. [ ]). Also, the 
choice of a local threshold criterion emphasizes the relation of these phenomena to other 
areas in nonlinear PDEs, where related phenomena occur, but the concept of black holes 
is absent. 

Critical behavior of the kind originally found by Choptuik is usually referred to as 
type II, because of its formal correspondence with type II phase transitions of statistical 
physics. A different type of critical solutions at the threshold of black hole formation, 
corresponding to type I phase transitions, is provided by unstable static configurations 

- like those found by Bartnik and McKinnon [BM88]. 

Linear perturbation calculations of such critical solutions revealed exactly one un- 
stable mode, which confirmed their interpretation as intermediate attractors in the lan- 
guage of dynamical systems. Critical phenomena in general relativity are reviewed in 
[Gun98, Biz96], including discussions in terms of phase transitions and renormalization 
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group techniques familiar from statistical physics. 

A massless scalar field in spherical symmetry exhibits type II critical collapse (there 
are no regular stationary or time-periodic solutions). Type II critical solutions have 
been found to exhibit continuous or discrete self-similarity in the past lightcone of the 
singularity. In our case, the critical solution is known to be discretely self-similar (DSS), 
and has been constructed directly as an eigenvalue problem [Gun95]. 

A spacetime is said to be DSS [Gun99] if it admits a discrete diffeomorphism Oa 
which leaves the metric invariant up to a constant scale factor: 



where A is a dimensionless real constant and n £ N. 

We choose scalar field critical collapse in spherical symmetry for several reasons: 
the model is very well studied and we can compare with a large amount of previous 
numerical and analytical results. Furthermore, the model is also very demanding: The 
value of the echoing period in the DSS critical solution is A ^ 3.44, which is quite larger 
compared to many other models. Note that larger values of A make it more difficult to 
resolve a large number of echos. 

Our numerical method is based on a characteristic initial value problem, i.e. we 
foliate spacetime by null cones. This allows for a very efficient evolution system and 
simplifies the study of the causal structure of the solutions. In spherical symmetry, 
caustics are restricted to the center of symmetry, so we do not have to deal with the 
dynamical appearance of caustics, which causes potential problems for characteristic 
initial value problems in higher dimensions. 

The numerical approach used in the compactified code[PHA05] is based on the 
"DICE" (Diamond Integral Characteristic Evolution) code, which has been documented 
in [HLP"'"00]. It mixes techniques from previous work of Garfinkle [Gar95] and the 
Pittsburgh group [GW92b, Win05], in particular we follow Garfinkle in moving along 
ingoing null geodesies to utilize gravitational focusing for increasing resolution in the 
region of large curvature. Furthermore, compactification methods are well studied and 
relatively straightforward to implement in characteristic codes. 

An important aspect of our compactified characteristic evolution scheme is that 
at late times our null slices asymptotically approach the event horizon, see Fig. 1.1. 
Essentially this is because our coordinates can not penetrate a dynamical horizon [ \K02, 
AK03, AK04, Hay94b] (they become singular at a marginally trapped surface, e.g. at an 
apparent horizon), which is spacelike if any matter or radiation falls through it and null 
otherwise [ E73, AK04]. Note that the dynamical horizon is contained inside of the 
event horizon, and the outermost dynamical horizon approaches (or coincides with) 
the event horizon at late times, assuming cosmic censorship holds. This fact makes 
our approach in some sense complementary to previous critical collapse studies, which 
were not adapted to the asymptotic regime. 

To further investigate physical quantities such as the mass function when a dynami- 
cal horizon forms, we employ an uncompactified double-null code which can penetrate 
dynamical horizons. This code is based on the work of Hamade and Stewart [ 1S96] 
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and improvements by Harada and Carr [HfirO'^]. The code has also been developed to 
analyze collapse problems in 2 + 1 dimensional gravity. 

In this thesis we focus on those aspects of critical collapse which are associated with 
global structure, and in particular the phenomenology seen by asymptotic observers. 
For the sake of completeness. We will also discuss other well-studied aspects such as 
mass scaling or local DSS behavior in the past self-similarity horizon. 



1.1 Organization of this Thesis 

This thesis is organized as follows: in chapter 2 we discuss characteristic surfaces of 
hyperbolic PDE's, mention different horizon concepts for black holes, and introduce our 
geometric setup, which is based on Bondi-type coordinates and double-null coordinates 
in spherical symmetry. In this context we also state gauge and regularity conditions at 
the origin. 

Chapter 3 we present evolution systems for the Einstein-massless scalar field prob- 
lem in both Bondi and double-null coordinates. We discuss a compactification scheme 
which introduces the Misner-Sharp mass-function as an independent evolution vari- 
able, which renders our evolution system regular at null infinity. We also mention 
physical quantities of interest, such as the Bondi-mass and Ricci scalar curvature. 

Our numerical algorithms are presented in chapter 4. In addition to the numerical 
schemes for the Bondi and double-null codes, we discuss issues of mesh-refinement, 
accuracy, convergence and the satisfaction of constraint equations in numerical evolu- 
tions. 

Chapter 5 deals with critical phenomena in gravitational collapse. We introduce the 
concept of continuous and discrete self-similarity, the dynamical systems picture and 
present a standard derivation of the mass-scaling law. 

In chapter 6, we introduce quasinormal modes and tails and present heuristic results 
for the presence of the least scalar damped QNM of a Schwarzschild black hole in 
critical collapse evolutions. We also study power-law tail exponents in different regions 
of numerical spacetime. 

In chapter 7 we give a detailed discussion of our numerical results pertaining to 
critical behavior, which include a comparison of the radiation signal at future null 
infinity with a heuristic estimate based on self-similar scaling and the accumulation of 
matter exterior to the past self-similarity horizon via backscattering. 

Our results and conclusions are summarized in chapter 8. 

In Appendix A and B we list tensor components in Bondi and double-null coor- 
dinates. Appendix C discusses numerical methods, starting from discretizations of 
ordinary differential equations and presenting a heuristic error analysis for the NSWE- 
algorithm. Finally, appendix D mentions the basic methodology which underlies con- 
vergence tests. 
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Figure 1.1: A Penrose diagram of a typical collapse spacetime. Shown is our numerical null 
grid which extends to future null infinity J^^. The grid consists of the null slices u - const 
and ingoing radial null geodesies v = const. Evolution slows down in the vicinity of the future 
event horizon We also indicate lines (a) r = const < 2Mf, (b) r = const = 2Mf and (c) 

r = const > 2Mf, where Mf is the final black hole mass. 
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1.2 Conventions 

We use the metric signature ( — h ++) and work in "geometrized units" G — c = 1. 
Spacetime indices are Latin. Space indices are denoted by Greek letters. Angular 
indices or indices of two-dimensional tensors are denoted by capital Latin letters. 
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Geometric Setup 



2.1 Some Notes on Hyperbolic PDEs 

In this thesis we will essentially be dealing with the numerical solution of nonlinear 
■partial differential equations (PDEs) on characteristic manifolds. Since this characteristic 
approach is entirely different from and not as common as the standard spacelike Cauchy 
problem, it is well worth the effort to investigate how characteristics and associated 
concepts feature in the theory of PDEs. 

A quasilinear (hyperbolic in our case) partial differential equation (PDE) of 2"** order 
in n variables Xi, x„ for an unknown function m(x,) can be written: (See Refs. [Garb4], 
[CH62]) 

" d^u " du 
L{u) = y ajk{du,u,Xi)- — h > bj{du,u,Xi)- \- c{du,u,Xi) + d - 0, (2.1.1) 

and thus defines the differential operator L{u). 

A submanifold of dimension n — 1 and equation (p{xj) - Ois said to be a characteristic 
manifold (surface) if its normals (pj fulfill the first order PDE^ 



j,k=i 



The Cauchy problem for a 2"'^ order hyperbolic PDE can only be well posed if the 
manifold on which the initial data is specified is not a characteristic manifold. If it 
were posed on a characteristic manifold, the PDE would only be first order in time, L(w) 
being an "interior" (i.e. acting purely in the characteristic surface) differential operator, 
and the 2"'' derivative transversal to the initial data surface could, in this case, not 
be determined. Thus, characteristic surfaces are exactly those, on which the Cauchy 



^The eikonal equation of geometric optics is of that form. 
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problem is not well-posed. Note that in the quasilinear case the characteristic condition 
depends on the initial data, since the aj^ are, in general, also functions of u(x,). 

Characteristic surfaces are, in turn, generated by characteristic rays or bicharacteristics, 
which satisfy the following equation 

)c=l 

It can be shown that the bicharacteristics of general relativity are null geodesies which 
generate null hypersurfaces. We will show in section 2.1.1 that the bicharacteristics 
of a massless scalar field which is minimally coupled to gravity are null geodesies of 
spacetime. Therefore, since gravitational disturbances, and, in our case, disturbances 
in the massless scalar field, are propagated along null geodesies, this makes a case 
for introducing a foliation of spacetime based on a family of non-intersecting null 
hypersurfaces. It is also worth noting that these disturbances need not be smooth in 
general, i.e. shock waves are possible in some matter models (e.g. perfect fluids). 
As a simple example consider the flat space wave equation 

Utt = UXX + Uyy + Uzz (2.1.4) 

The characteristics then have to fulfill 

<^] = cg^c^'y + cg. (2.1.5) 

A special solution is the characteristic cone 

(/) = (f - t^f - (X - x^f - (y - yo)' - (z - zo)' = 0, (2.1.6) 

which represents spherical wavefronts propagating from a source located at (xq, i/O/Zo). 
The generators of the cone, which are bicharacteristics, may be represented as light rays 
perpendicular to the wavefronts. 

It is illuminating to compare the spacelike slices of the 3-1-1 formalism with a char- 
acteristic foliation based on outgoing null hypersurfaces for the scalar wave equation 
□(/) = in Minkowski spacetime with polar spherical spatial coordinates (r, Q, (p). 

In 3 -I- 1 the result is well known: 

11 1 

dtt(p = -^dr{r%(p) + . deisin Odgcf) + -^^d^^(p. (2.1.7) 
r^ ^ ^ r^smO ^ ' > r^^ sm 9 ' ^ 

This equation is second order in time, i.e. both the field and its time derivative must be 
specified on an initial slice in order to uniquely determine the future evolution. 

In contrast, let us introduce a coordinate system {u,r,6,(p) based on null cones 
u = const, which are generated by a two-dimensional set of null rays, labeled by the 
angles 6, (p and with r a radial parameter along the rays. These coordinates will be 
discussed in detail in section 2.3.2. The line element is 

ds^ = -du^ - Idudr + r^dCl^, (2.1.8) 
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so that the wave equation then becomes 
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2{r(P),ur = {r(t)),rr ' ji^\r(\>) (2.1.9) 

where is the angular momentum operator. This equation is only first order in the 
retarded time u, so only the value of the field needs to be specified on an initial slice in 
order to uniquely determine the future evolution. However, the domain of dependence 
only extends to the inside of the cone; i.e. the initial null surface is not a Cauchy surface 
for the spacetime. 

2.1.1 The Bicharacteristics of a Scalar Field Coupled to Gravity 

Let us consider the matter field equation for the massless scalar field minimally coupled 
to gravity in spherical symmetry and assume that we can foliate spacetime by null hy- 
persurfaces which are labelled by a coordinate u. We will show that the bicharacteristics 
of the scalar field are null geodesies of spacetime (see [d'I92b]). 
The matter field equation is of the form of equation (2.1.1) 

Ugcp = VV^c/) = ^^dad^cp - ^""rjccp, (2.1.10) 

and its principal part (the terms which contain the highest derivatives) is 

P[Ug(p\=^''dad^(P. (2.1.11) 

The coeffient matrix which appears here, the Uj^ of (2.1.1), is just the inverse spacetime 
metric in the (u,r)-submanifold. According to equation (2.1.2) the u - const hypersur- 
faces are characteristic manifolds since they fulfill 

/^V«MVbM = 0. (2.1.12) 

By definition (2.1.3) the bicharacteristics which generate the u - const surfaces satisfy 
(with A being a parameter along the bicharacteristics) 

dA 

i.e., they are the orbits of the vector field Vw. 

Now we will prove that these bicharacteristics are, in fact, null geodesies. We take 
the covariant derivative in the direction of the vector field Vw of the last equation and 
find 

dx" \ ^^ij 



^=/^V,M, (2.1.13) 



Thanks to the symmetry of the connection this last expression is equal to 

/VwVfcV.M = i/^Vfc (V'uVcu) = 0, (2.1.15) 

where we have used that Vw is null. Thus, the bicharacteristics of the massless scalar 
field turn out to be null geodesies. This intuitively makes sense, since the scalar field 
propagates at the speed of light. We have also shown that Vu is an affine null vector, 
that is, A (chosen by equation (2.1.13)), is an affine parameter along these null rays. 
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2.2 General Properties of Black Holes 



In the following we give some textbook definitions of a black hole, both from a global 
and a local point of view [Poi04, Wal84]. 

The causal past of an event p, denoted by /~(p) is the set of all events that can be 
reached from p by past-directed curves, either timelike or null. The causal past of a set 
of events S, /~(S) is the union of the causal pasts of all events p € S. 

The black hole region B of the spacetime manifold M is the set of all events that do not 
belong to the causal past of future null infinity: 

B^M-J-{J^*). (2.2.1) 

The event horizon EH is defined to be the boundary of the black hole region: 

EH^dB^d iriJ^^)) . (2.2.2) 

Because the event horizon is a causal boundary, it must be a null hypersurface generated 
by null geodesies. Note that this definition of a black hole depends crucially on the 
notion of future null infinity which can be defined for asymptotically flat spacetimes. 
Moreover, to apply this definition one has to know the entire future development of 
the spacetime in question. In practice, it is often more convenient to use a quasi-local 
definition of a horizon. 

A trapped surface is a closed, spacelike two-dimensional surface S such that for both 
congruences (ingoing and outgoing) of future-directed null geodesies orthogonal to S, 
the expansion 6 is negative ever}rwhere on S. 

The trapping ?zorzzon[iiay94a] is the three-dimensional boundary of the region of 
spacetime that contains trapped surfaces. 

The two-dimensional intersection of the trapping horizon with a spacelike hypersur- 
face E (Cauchy surface) is called an apparent horizon. The apparent horizon is therefore 
a marginally trapped surface - a closed two-surface S on which one of the congruences of 
future-directed null geodesies orthogonal to S has zero expansion. ^ Note that the inter- 
section of an outgoing null hypersurface with the trapping horizon can be an apparent 
horizon, or a three-dimensional subset of the trapping horizon. 

The trapped region ^ of a spacelike h}^ersurface E is the portion of E that contains 
trapped surfaces. 

Under certain technical conditions the black hole region contains trapped surfaces. 
While for stationary black holes the trapping horizon typically coincides with the event 
horizon, it lies within the black hole region in dynamical situations (unless the null 
energy condition is violated.) We will not distinguish between the two-dimensional 

trapping horizon is the closure /? of a three-surface H foliated by marginal surfaces on which 0_|h 7^ 
and # 0. a marginal surface is a spatial two-surface S on which one null expansion vanishes, fixed 

as 0+ls = 0. The trapping horizon and marginal surfaces are said to be outer if X-0+|h < and future if 
0_|h < 0. The concept of dynamical horizons[ ] is closely related to trapping horizons. 

^For technical details see [Wal84] theorem 12.2.5. 

*For details see [ ] proposition 12.2.3. 
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apparent horizon and the three-dimensional trapping horizon and will refer to both as 
the apparent horizon. We also use the term "trapped region" in a different sense as 
defined above, namely we refer to a region of spacetime for which each point lies on a 
trapped surface [Hay94a]. 

In the following we will restrict ourselves to future outer trapping horizons (FOTH) 
and future marginally outer trapped surfaces (MOTS), which require 0+ = and 6- < 0, 
where + and - refers to the outgoing and ingoing null congruences, respectively. A 
FOTH is either spacelike or null, being null if and only if the shear of the outer null 
normal as well as the matter flux across the horizon vanishes [ * ' ']. 

The apparent horizon is the outermost (future) marginally outer trapped surface 
(MOTS) in a spacelike hypersurface E. If one looks for globally outermost MOTS in a 
Cauchy evolution, they can "jump" [AMS05]. This does not happen for an outgoing 
null foliation. 

2.3 Bondi Coordinates in Spherical Symmetry 
2,3.1 Spherical Symmetry 

By definition, the metric of a spherically symmetric spacetime possesses three spacelike 
Killing vector fields, which form the Lie-algebra of SO(3), the rotation group in three 
dimensions. The Killing vectors fulfill the following commutation relations 

[v„Vj]^e1.V,, (2.3.1) 

where the e^. denote the structure constants of SO(3). The commutators close (that is, 
the commutator of any two fields in the set is a linear combination of other fields in 
the set), so that by Frobenius' theorem (see [VVal8 ]) the integral curves of these vector 
fields "mesh" together to form submanifolds of the manifold on which they are defined. 
The orbits of SO(3) are spacelike 2-spheres and the spacetime metric induces a metric 
on each orbit 2-sphere. 

The general form of a spherically symmetric metric is (see appendix B of [HE73]) 

ds^ = dz^it, r) + Y^{t, r)dD^iO, cp), (2.3.2) 

where dz'^ is a two-dimensional metric of indefinite signature ( — h) and dQ.^ — dQ^ + 
sin^ Odcp^ is the canonical metric on S^. The coordinates (f, r, 0, (p) are chosen so that the 
group orbits of SO(3) are the surfaces t,r — const and the orthogonal surfaces are given 
hy 6,(p - const. 

23.2 Bondi Coordinates 

In the following, we define Bondi coordinates and derive the form of the spacetime 
metric, (see [dT92b] and [Pap74]). We start from the general form of a spherically 
symmetric metric given in equation (2.3.2). We chose 

x° = M = const (2.3.3) 



Bondi Coordinates in Spherical Symmetry 



11 



to be a family of non-intersecting, outgoing null hypersurfaces, i.e. 
g~^{du, du) — 0. We know from section 2.1.1 that these null hypersurfaces are generated 
by null geodesies (the bicharacteristics of the surfaces). This gives us the possibility to 
define a second coordinate in a purely geometrical way. We choose 

= r (2.3.4) 

to be a radial parameter along the outgoing null geodesies that foliate the u = const 
surfaces, i.e. on each such null ray we have u = const and and - const. 

The remaining coordinates are chosen to be polar angles, coinciding with the usual 
flat space definitions for 

x^ = e,x^^(p (2.3.5) 

near the center of spherical symmetry, as we will later choose u to be central proper 
time. This completes our coordinate system to (w, r, 6, cp). 

Now, we derive the form of the spacetime metric in these coordinates. An outgoing 
null ray (light ray) is a coordinate curve 

u = uq, Q - Qq, (p = (po all constant (2.3.6) 

where r is varying, or 

du^de^dcp^O (2.3.7) 
The tangent vector to such a curve x"(r), 

d Y dx' 



must be parallel to the normal vector to the u = const hypersurfaces 

Vm = /''VfoM = (2.3.9) 

since, in the case of a null hypersurface, the normal vector is also the tangent and the 
null rays generate the u = const hypersurfaces. Thus, we find 

g"" cc 6^, (2.3.10) 

which entails 

= 0, (2.3.11) 

or equivalently 

grr = 0. (2.3.12) 

We still have freedom in the choice of parametrization along the null rays. Following 
Bondi, we choose r to be an areal radial coordinate or luminosity distance, so that the 
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surface area of r = const, u = const 2-spheres equals Anr^. ^ This choice can be enforced 
by the condition 

sin^ 9. 



(2.3.14) 



Note that in spherical symmetry, the areal r-coordinate is defined in a geometrically 
unique way ([iiay96]). It is also possible to compactify the radial coordinate; we 
discuss this in detail in section 3.6. 

Putting these requirements together, we find that the line element reduces to the 
form 

ds^ = guu{u, r)du^ + Iguriu, r)dudr + r^dQ^. (2.3.15) 
Since (^|"," |"^' ) is Lorentzian, we can choose to be timelike, i.e. 

d d 



du ' du 



<0. 



(2.3.16) 



We write the metric in the following way, adhering to Bondi's conventions for the 
names of the metric functions 



ds^ = -e^Piu,r)Y^au2 _ 2e^PM^udr + r^dn\ 



(2.3.17) 



where V > is assumed in order that ^ be timelike. 



So the metric and inverse metric are 



-e^/ 
n a 








sin^ i 



(2.3.18) 



^Other possible choices include an affine parametrization of the outgoing null rays, as used by Newman 
and Penrose in their work on gravitational radiation [NP62]. An affine radial coordinate entails by 
definition 

which leads to 

= drd";. + r,.^b'i = Ti.,.. 

For null coordinates in spherical symmetry the only nonzero contribution comes from 

K=g"'gur,r. 

A vanishing of g'"' and thus diverging of g,,,- would make the metric singular, therefore the condition to be 
satisfied is 

gur,r = 0. 

According to the general form of a spherically symmetric metric given in equation (2.3.2), the line element 
must then be of the following form 



ds^ = -guu(u, r)du^ - Idudr + f{u, r)dD} 



(2.3.13) 
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and 



-e-^P 

-e-2p e-^^Y. 

' ;3 

'o 









(2.3.19) 



r2 sin^ Q 

respectively. 

The sign of gur is due to our use of a retarded null-coordinate w. In section 2.3.4 we 
will choose u to be proper time at the origin which implies that the metric reduces to 
the following flat metric at the origin 

^^)\at = ~ 2dwdr + ?dQ^. (2.3.20) 

One can also check that this is consistent with the signature ( — h ++) by introducing a 
time coordinate i = w + r in flat space, so that 



(2.3.21) 



For purposes of comparison, we can also introduce an advanced null-coordinate v = t+r 
in flat space, and find 

dsj^^^ = -dv^ + Idvdr + r^dCl^. (2.3.22) 




Figure 2.1: Bondi coordinates (m, r, d,(p). u is the proper time at the center of spherical symmetry, 
r is an areal radial coordinate. A regular center is assumed at r = 0. 

In the context of Bondi coordinates, we often write partial derivatives of some 
function f{u, r) with respect to u and r as / and /', repectively. 
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2.3.3 Physical Interpretation of the Bondi Metric 



To improve physical understanding, it is worthwhile to analyze the Bondi metric in 
some further detail. 

First, let us take a look at the tangent and normal vectors as shown in table 2.1. By 
construction, we have 

= -e^N'u (2.3.23) 

i.e., (^^j lies in the u=const null surfaces. Or, to put it differently, these null rays 
generate the u — const null hypersurfaces. 

Observe also that if V is positive, then r = const is a timelike surface and (d/du)" is 
also timelike. lfV-0, both are null, and if y < both are spacelike. Correspondingly, 
V^r is spacelike, null or timelike. 





(1,0,0,0) 


Vm = 


(0,-e-2^0,0) 


V^mVm = 


null 




(0,1,0,0) 


V«r = ( 


_e-2/5^e-2/^Z 0,0) 


V^rVr = > 


spacelike 


(^) 


= (1,0,0,0) 


(^) - 

V du )a 






timelike 


(ir 


= (0,1,0,0) 




(-e2^ 0,0,0) 


(ir(ii='' 


null 



Table 2.1: Overview of tangent and normal vectors to coordinate hypersurfaces and coordinate 
lines. Here we rely on the assumption that V > 0. 



The particular choice of the algebraic form of the metric components serves the 
simplicity of the resulting Einstein equations. The functions |S and V also have sim- 
ple physical interpretations. As will be shown in section 2.3.4 |S measures the redshift 
between an asymptotic coordinate frame of reference and the center of spherical sym- 
metry. V is the analog of the Newtonian potential and contains the Bondi-mass in its 
asymptotic expansion. 

We define ingoing and outgoing null directions /+ by 

= (l,-^y,0,0), (2.3.24) 

and 

II = - V'm = (0, e-^^, 0, 0), (2.3.25) 

which is an affine null vector, i.e. l^yj^. = 0. For an affine null vector the null expansion 
@ can be defined by (see chapter 9 of [Wal84]) 

e = yj". (2.3.26) 

We find 

e+ = (2.3.27) 
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Figure 2.2: This figure shows how the u = const null slices approach an apparent horizon in the 
late part of a near-critical time evolution. Observe how the first slices shown are still close to a 
flat space null-cone, whereas the last slices extend almost in a vertical direction, i.e. 



dv 
dr 



If ©+ > along an outgoing null ray, the surface area of 2-spheres increases and 
neighboring null rays diverge. Since 0+ is positive definite unless |S diverges, the 
occurrence of marginally outer trapped surfaces, 0+ = for some sphere, cannot happen 
in outgoing Bondi coordinates. Thus, we cannot penetrate apparent horizons (see 2.2). 

Nevertheless, as we come close to an apparent horizon, our slices will asymptote to 
it (see figures 2.2 and 2.3) and ultimately, as 0+ ^ 0, r will cease to be a good coordinate; 
it will no longer be a monotonically increasing function of the advanced null coordinate 
V which is defined by ^ = /_, i.e. the curves v - const are the ingoing nuU-geodesics 
and V = r on the initial slice. (In flat space v = t + r.) Tiny numerical errors may lead to 
the bending over of the slices and thus crash the code. 

The radial ingoing null geodesies are given by integration of the equation 



±r(u) = -YMm 
du ^ ' 2r{u) ' 



(2.3.28) 



Accordingly, the total derivative of any scalar quantity / with respect to u along an 
ingoing null geodesic is given by 



(2.3.29) 



Since this is a Lie-derivative with respect to a non-coordinate vector field, it will in 
general not commute with partial derivatives like J^. 



2.3.4 Gauge and Regularity Conditions at the Origin 

Assume that there is a regular center at r = (no eternal black or white holes or naked 
singularities), so that we can designate a finite geodesic line segment of a central observer 
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0.2 0.4 0.6 0.8 1 

l+r- 

Figure 2.3: This figure shows the same situation as depicted in 2.2, but for a near-critical 
compactified run. The initial slice is by definition a straight line. All slices extend to J'^ at x = 1. 
Late slices exhibit the formation of an apparent horizon close to the origin. 



at r = as our spatial origin. We fix the gauge such that the family of outgoing null 
cones emanating from the center is parametrized by the proper time ^, u, of this central 
observer. This choice is taken since we are interested in the "echoing region" near the 
center of spherical symmetry, where oscillations in the fields will pile up if we are near 
to the "critical solution". Another common gauge is to choose u as the proper time of 
an observer at future null infinity, J^"*". This asymptotically flat time coordinate is the 
usual choice if one is interested in asymptotics. 

Choosing u as proper time at the origin implies the condition 

e^^- = 1, (2.3.30) 

r 

while a regular center, as defined by [Hay96], requires g~^{dr,dr) - 1 = 0{r^), which in 
turn entails 

e-^f^Y. ^l+0{r^). (2.3.31) 
Together theses requirements imply local flatness and regularity of the metric, that is 

e'P{uo,r)^l+0{f^) 

V , (2.3.32) 

-(uo,r) = l+0(r2), 

at a fixed retarded time Uq. These smoothness conditions are rederived in section 2.3.4. 



^ With our choice of Lorentz signature (- + ++), the proper time of a timelike curve C with tangent T" 
and curve parameter t is given by 

T 



J {-gnbT'T')'" dt. 
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Relation between Central and Bondi Time 



In the following, Uq will always denote central time - proper time at the center of spheri- 
cal symmetry, whereas we denote the asymptotically flat Bondi time by Mb- The relation 
between the two time coordinates can easily be derived by comparing with an asymp- 
totically flat (standard Bondi) frame (see section 3.6.1). For large r the line elements 
behave as 



ds^ = - (e^^ + 0{r-^)) duc [(e^^ + Oir'^)) duc + 2dr\ + r^dO} 
ds^ - -duJ' - Iduudr + r^dQ^ 



which yields in the limit r — > oo 



du^ - e™duc, 



(2.3.33) 
(2.3.34) 

(2.3.35) 



where H - |S(Mc, oo) and lim^^oo j{uc,r) = e^. 

The following more geometrical argument yields the same result: In general Mb = 
Wb(mc/'')/ but under the following assumptions we can show that = Mb(wc) only. We 
demand that 

<o. 



and 



VmbV^Mb = 0, 



(2.3.36) 
(2.3.37) 



i.e. Mb = const are null hypersurfaces. 
We find 



= V"uj,{uc,r)VaU^iuc,r) 



-21} 



IduAlVfdu, 



du 
\ dr 



r \ dr 



-2 



dtic 



which yields 

Thus Wb = Mb(mc) only. Now we have 

'A 



dr 



0. 



= /("c) 



Since is proper time at infinity. 



d_ d_ 
dtiu' dih 



plugging (2.3.40) in to the metric yields 

\duc' duc 



ducj 



-I, 



r 



2 AH 



(2.3.38) 



(2.3.39) 



(2.3.40) 



(2.3.41) 



(2.3.42) 
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Therefore, f = e and we have found the desired relation between central and Bondi 
time: 

= e™duc. (2.3.43) 

Redshift 

Consider a freely falling observer located at r = whose worldline is tangent to the 
timelike geodesic parametrized by proper time at the origin Uq with tangent (d/duc)". 
He sends out a light signal which propagates along the affinely parametrized outgoing 
radial null geodesic with tangent = —Vuc (As shown in section 2.1.1, k" is an affine 
null vector). The signal is received by an observer whose worldline is tangent to a 
timelike geodesic at r = oo, parametrized by proper time at infinity Wg with tangent 

id/du,T. 

The frequency (see Ref. [WaI84]) of an emitted wave is given by minus the inner 
product of the wave vector (which is just the rate of change of the phase of the wave) 
with the 4-velocity of the observer. The frequency of emission at the center is given by 



r=0 \o'"c 



= 1. (2.3.44) 

r=0 



The frequency of reception at infinity is, using the relation between central and Bondi 
time (2.3.35), 



(Ob = -/Cfl I ^ 



The redshift-factor is defined as 



a 

2H 



= e-^. (2.3.45) 



z = ^-l = e2«-l. (2.3.46) 

COB 

For small H we thus have z = 2H + 0{H^). 

Now, because of the hypersurface equations (see section 3.3) and gauge choice for 
u, the metric functions, jS and V, are both monotonically increasing in the direction of 
increasing r and both are positive definite. Thus H = jS(m, oo) > and 

COB < coc (2.3.47) 

always, which is just the condition for signals to be redshifted. In the case of horizon 
formation 0+ — > 0, however, H will tend to infinity, so that light rays can no longer escape 
to infinity and the redshift will diverge exponentially until our coordinates break down. 

A crude argument for the redshift is also apparent just from the relation (2.3.35). If 
H — > cx), a finite amount of central time Uq corresponds to an infinite amount of Bondi 
time Mb arid thus light signals originating from the center are infinitely redshifted. 
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Derivation of Regularity Conditions at the Origin of Spherical Symmetry 

In this section we apply the general procedure for deriving regularity conditions set 
forth by Bardeen and Piran in [ ]. Although the latter deals with axisymmetric 
systems, it is nevertheless worthwhile to study this approach in the simpler setting of 
spherical symmetry. 

A tensorial quantity is called regular at r = if and only if its components relative 
to Cartesian coordinates can be expanded in non-negative integer powers ^ of x, y and 
z; i.e. one demands the existence of a Taylor series expansion in Cartesian coordinates. 
We would like to derive what restrictions this assumption of regularity at the origin 
implies for the metric functions |S and V as functions of r. 

We will work in a regular {t, x, y, z) coordinate system with t defined as t - u + r. ^ 
Since we have taken ii to be proper time at the origin, t also measures proper time at 
r = 0. We also enforce local flatness at the center, i.e. the metric goes to the Minkowski 
metric as r tends to 0, 

lim^flfc = T]ah- (2.3.48) 

We will have to investigate which derivatives of the metric functions are allowed to be 
non-zero and which are zero by regularity. 

The relations between the two coordinate systems are: 

t = u + r 
X = r sin 6 cos (p 



y = r sin 6 sin ^ 
z = r cos 6 



u = t - r 

r - yjx^ + y^ +z'^ 



6 = arctan 



y 

= arctan - 



x^ + 1/2 



The Jacobian for this coordinate transformation is 

o'(M,r,g,(/)) 
d{t,x, y,z) 



First we need to calculate the components of the metric in {t, x, y, z) coordinates. 
Some interesting components as functions of the polar coordinates are 

gn = -e'^y (2.3.50) 
^ Otherwise the tensorial quantity or its derivatives will blow up at r = 0. 

^Note that the coordinates {u,x, y,z) are not regular: The vector fields ^, ^ and ^ all have a kink at 
the origin. 



1 
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_l 





r 
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r 
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y 


X 




x^+f 



■V 



_z 

r 

z 
r 



x^+y^ 



(2.3.49) 
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gtx = — sin 6 cos (p - e^^ sin 6 cos (p 



(2.3.51) 



g^^ = -e^^ Y- sin^ e cos^ (p + 26^^ sin^ 6 cos^ c^) + cos^ 6 cos^ + sin^ (2.3.52) 



gyy = -e^^ Y- sin^ e sin^ + 2e^^ sin^ 6 sin^ + cos^ 6 sin^ + cos^ (j) (2.3.53) 

= -e^P Y- cos^ + le^f^ cos^ + sin^ 6 (2.3.54) 

We will also need derivatives of the metric. Since the derivatives of gxx and gyy are 
particularly messy and we do not need to consider them, we leave them out. 



Stt^ 



V y sin cos (i) sin cos (6 
= -7j^^^,r- sin cos (^) + - ^^V,r (2.3.55) 



gtx,x 



sm fcos 



+ e 



2V 1 V, 

+ - + 



sin2 0cos2(^). (2.3.56) 



g^^^^ = ^-2j3,,e^^ - + 2^^ ^ + 4jS,^e2^ - — j cos^ sin cos c^) 

^ ( r ~ ~ ^)) ^ '^os (^) - ^ sin^ cos (p. (2.3., 

We make the following series ansatz for the metric functions jS and j: 



p = a + br + ct^ + 0{r^) (2.3.58) 

j=d + er + fr^ + 0{r^) (2.3.59) 

This yields 

^,r = b + Icr + 0{r^) (2.3.60) 

V,r = d + ler + 3fr^ + 0{r^) (2.3.61) 

e^^ = (1 + 2fl + 2fl2) + (lb + Aab)r + (2c + 2^^ + 4ac)r^ + 0{r^). (2.3.62) 
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We now impose the local flatness gauge condition (2.3.48). 



yields 

Similarly, we have 



lim^tt = -1 (2.3.63) 

)-->0 



£l(l + 2fl + 2fl2) = -1. (2.3.64) 



lim = lim e^^ sin cos (f) f- - l) = 0, (2.3.65) 

from which we find 

d = l, (2.3.66) 

and thus, it follows from (2.3.64) that 

fl = 0. (2.3.67) 

One does not need to consider any other equations, since the two constant terms in 
(2.3.58) have now been determined. Indeed, one can also check that lim,-^o gxx = 1 arid 
lim^^o gzz = 1 are fulfilled. 

The series that respect this gauge choice are 



.2 ^ ^(^3) 

l+er+fr^ + 0(r^) 
r 

J2.\ 



V 



^,r = b + let + 0{r^) (2.3.68) 
y , = 1 + 2er + 3fr^ + 0{r^) 
e^l^ = 1+ 2br + (2c + 2b^)r^ + 0{r^). 

We are now in a position to investigate limits (r — > 0) of the derivatives of the metric. 
The values should remain finite and be direction independent, i.e. independent of the 
spherical angles 6 and (p, since these are indefinite at the origin. This will then provide 
us with conditions on the series expansions of j3 and j. Consequently, we find from 
(2.3.55), the equation for gtt,x, that 

-2^,re'P^+e'^l^-MiO{r). (2.3.69) 

Applying the expansions (2.3.68), we find 



(l + 2br + 0{r^)) 
and eventually 



1 1 

- + e+ fr 2e - 3fr - 2& - 4cr - Iber + 0(r^) 

r r 



= 0{r) (2.3.70) 
e + 2b = 0. (2.3.71) 
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Similarly, we find from (2.3.57), the equation for gzz,x, that 

(-2/^,1 + 4/^, + f - ^) cos^ - J^) + ' (1 - sin^ 0) ^ Oir). (2.3.72) 



Here, we can set all terms proportional to cos^ 6 to zero and insert the series expansions 
(2.3.68). We find 

e + 2fr + Iber + 0{r^) = 0{r), (2.3.73) 

which leads to the condition 

e = 0. (2.3.74) 
Combining the two conditions (2.3.71) and (2.3.74) we recover 

b = (2.3.75) 
e = 0. (2.3.76) 

We did not consider gtx,x, since it is already well behaved because of the gauge 
condition alone. For gxx,x the calculation would have been more involved but would 
have led to the same result. 

Summarizing, we have found that the metric functions (5 and j guarantee that the 
metric is regular at the center of spherical symmetry if they satisfy, at a fixed retarded 
time uq, 









= 1 + 0{r^) 



(2.3.77) 



These regularity conditions are consistent with the hypersurface equation (3.3.3), V = 
e^f^, derived in chapter 3.3. 

The conditions we have derived are sufficiently accurate for our purposes, though 
one could consider higher derivatives of the metric. We will see from the following 
simple argument that, in fact, all odd derivatives of the metric functions have to vanish 
at the origin. 

Without loss of generality, define linear distance I by 

i-r if z<0, , , 

l = z = \ 2.3.78 
1+r if z>0. 

We demand that the metric does not have a kink as the origin is crossed along a ray 
parametrized by linear distance /. If a metric function /(r) contains a linear term in its 
Taylor expansion around the origin, then 
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where C is a constant. 



1 



-C + 0(r) if r < 
fir) =\ ^ ^ 2.3.80) 

\+C + 0{r) if r>0 ^ ' 



has a jump discontinuity and is thus not regular. It will be regular provided that the 
constant C vanishes. This is equivalent to the condition that 



= (2.3.81) 

r=0 



Note that this argument can be extended to higher derivatives. Roughly speaking, each 
time one applies ^ to a function f{r), f{r) reverses its sign for r < 0. Therefore, all 
derivatives of an odd order have to vanish at the origin, so that a regular function of r 
will have a Taylor series which contains only even terms in r: 

fir) = c + C2r^ + C4r^ + .... (2.3.82) 

In addition, we can also derive the Taylor series expansion for the scalar field by 

making use of the hypersurface equation (3.3.2) = 2nr {<p,i-j and find that there is no 
restriction on cp: 

(Pir) = cPo + (pir + Oir^). (2.3.83) 

In this respect it is also interesting to view this last fact from a different perspective: 
When working in polar coordinate systems, the d'Alembertian usually contains terms 

similar to ^dr- Regularity of the field (p at the origin then amounts to requiring that ^ lr=o 
vanishes. But using a null coordinate u makes a difference. There is now a rivaling 
-jdu term in the d'Alembertian (A.0.22) which happens to cancel the other term at the 
origin. 



2.3.5 Historical Notes on the Bondi Coordinate System 

The original motivation for Bondi and his co-workers (see [HB62]) in developing this co- 
ordinate system was to analyze radiation from an isolated system. To do so, one needed 
to deal with expansions in negative powers of a radial distance for various quantities. 
Such a treatment is, however, made impossible by the appearance of logarithmic terms 
in r, as is the case, in Schwarzschild coordinates. 

Later, the conformal compactification of null infinity allowed a geometrized 
description of the asymptotic physical properties of radiative spacetimes. This novel 
technique also proved very beneficial for the study of gravitational radiation, see e.g. 
[Fra04]. 

One demands that in a region from J^"*" inwards the coordinate system should be 
non-singular. Farther in, null rays will, in general, focus and cross, forming caustics 
and causing the coordinates to break down and become singular. ^ 

'Except for the case where there exists a secor\d regular center, where rays may caustic, this cannot 
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The Bondi-like coordinates defined in section 2.3.2 cover all of spacetime provided 
there are no caustics or singularities and J'^ exists. 



1A Double-Null Coordinates 

As we have seen in section 2.3.3, outgoing Bondi coordinates cannot penetrate apparent 
horizons, because the areal radial coordinate r fails to be a good coordinate when an AH 
forms: Vr — > 0. A coordinate system which does not suffer from this deficiency and is 
suitable for spherically symmetric critical collapse evolutions are double-null coordinates. 

In addition to the retarded null-coordinate u which is constant on outgoing light 
rays, we now introduce an advanced null-coordinate v, which is constant on ingoing 
light rays and replaces r as a coordinate. Along with the standard polar coordinates on 
the group orbits of SO(3), S^, the coordinate chart becomes (w, v, 6, cp). The line element 
is then of the following form 



ds^ = -a^{u,v)dudv + r^{u,v)dCf-. 



(2.4.1) 



The area radius r is now a metric function that depends on the null coordinates. The 
double-null form is natural in the sense that each two-sphere possesses two preferred 
normal directions, the null directions d/du and d/d^. 

The components of the metric and inverse metric are given by 



and 
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and the metric determinant is 



detg = - 



a'^r^ sin^ i 



(2.4.2) 



(2.4.3) 



(2.4.4) 



happen in spherical symmetry, because the lightrays are 'running' after each other but never touch. A 
second regular center is however ruled out by the assumed existence of future null infinity .y"*". The 
only caustic that is present in spherical symmetry is a point caustic at the center itself, which has to be 
handled by regularity conditions on the metric and the matter fields. Of course, for nontrivial topologies 
the situation is entirely different: e.g. on the Einstein cylinder lightrays may form caustics. 

^"For the existence of which is a part of conformal infinity, the spacetime metric has to satisfy certain 
falloff conditions, i.e. it must tend to a flat metric sufficiently fast as one goes to infinity in an outgoing 
null direction. In a coordinate independent formulation a spacetime is said to be asymptotically flat, if the 
physical spacetime can be mapped into a new, "unphysical" spacetime via a conformal isometry which 
has to satisfy certain, very technical conditions. Details can be found in chapter 11 of Ref. [Wal84]. 
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A relation between Bondi coordinates and double-null coordinates can be estab- 
lished by a coordinate transformation (see Ref. [ ]) of the form 

'V 



dv = e 



2a 



(^jdu + 2drj, (2.4.5) 



with an integrating factor a{u, r). 
The line-element is 



ds^ = -e^^ {^-du + 2drj du + r^dCl^ 
= -e^^-"dudv + r^dQ^. 



(2.4.6) 



In the new coordinates r - r{u, v) and we can introduce a metric function a{u, v), so that 
a^{u,v) = e2/5(".'(",^'))-«(",'-(",^')) and the line-element becomes equation (2.4.1). 
Since dr = {dr/du)\vdu + {dr/dv)\udv we have the relations 



du 
dv 



V 

\e-^" (2.4.8) 



In the context of double-null coordinates we denote partial derivatives of some 
function f{u, v) with respect to u and vhy f and respectively. 

2.4.1 Null Expansions 

To compute marginally trapped surfaces we compute the null expansions (see Ref. 
[Wal84]) 

0± = V„/l, (2.4.9) 

where the affine ingoing and outgoing null vectors Z+ are defined through their covectors 
as 

ih), := -Vau (2.4.10) 

(L)„ := -Vav. (2.4.11) 

This yields 

/« = = (0, 2/fl^ 0, 0)" (2.4.12) 

/« = -g^-^ = (2/fl^ 0, 0, 0)" (2.4. 13) 

It is easy to check that the /+ are indeed affine null vectors, that is /^V^/^ = and that 
the null expansions 6+ are given by 

4r' 

0+ = — (2.4.14) 

4r 

0- = — . (2.4.15) 
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The signs of the null expansions are geometrical invariants, while their actual values 
are not [ ]. For future marginally outer trapped surfaces (MOTS) (see section 2.2) 
we need to have G+ — while G- < 0. Since the denominators in equations (2.4.14) and 
(2.4.15) are non-negative, the sign of the null expansions is determined by the quantities 
r' and r and it suffices to check the conditions r' = while r < 0. 

The possible presence of multiple MOTS on an evolution hypersurface depends 
on the slicing. For outgoing null slices, there cannot be more than one MOTS on 
a slice, whereas for spacelike slices there can be. The geometric property of future 
outer trapping horizons being spacelike or null (the latter only if there is no more 
infalling matter when following the null geodesic generators further outwards), but 
never timelike, applied to the intersection of an MOTS with an outgoing null-slice 
shows that we cannot have multiple MOTS on such a null-slice, which is equivalent to 
the function r'{u = const, v) having at most one zero. 

2.4.2 Gauge Choice 

The form of the metric is preserved under diffeomorphisms u — > U{u) and v — > V{v). 
We fix this residual gauge freedom present in the null coordinates u and v by choosing 

• r = to be at M = 

• r{u = 0, t;) = v/2. 
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Chapter ^ 

The Continuum Problem 



We consider the coupled spherically symmetric Einstein-massless scalar field system, 
both in Bondi and in double-null coordinates. After deriving the Einstein equations 
and the curved space wave equation for the scalar field, we analyze the evolution 
systems and formulate the characteristic initial value problem to be solved. We also 
discuss compactification for Bondi coordinates and highlight some quantities of physical 
interest which we will use in the numerical evolution codes described in the following 
chapter. 

3.1 Einstein's Equations 

To derive Einstein's equations we start with the action functional which consists of the 
Einstein-Hilbert action plus the matter action: 

where R denotes the Ricci scalar and the Lagrangian density of the massless scalar 
(Klein-Gordon) field is given by 

^ = -^/''V«(/)V,(/). (3.1.2) 

Variation of the total action (3.1.1) with respect to the spacetime metric yields the 
Einstein equations, 

Gah = %nTah (3.1.3) 

where 

Tab = V«(/)Vfc(/) - ^^«feV,(/)V> (3.1.4) 

is the energy momentum tensor of the massless scalar field obtained via variation of the 
matter action with respect to the spacetime metric. 
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To complete the system, variation of the matter action with respect to the scalar field 
(p yields the matter field equation, the curved wave equation for the massless scalar 
field 

Dg(p = VVaCp = 0. (3.1.5) 

The Einstein and Ricci tensors in outgoing Bondi coordinates (w, r, 6, (p) are given 
in appendix A along with the Christoffel symbols, while the tensor components for 
double-null coordinates are given in appendix B. 

3.2 Hierarchy of Einstein Equations in Bondi Coordinates 

In spherical symmetry, the four algebraically independent Einstein equations, the (u, u), 
(u, r), (r, r), and (0, Q) components of equation (3.1.3), are not differentially independent. 

Essentially, we have four equations for two unknown metric functions, |S and V. In 
principle, one could freely choose two of these equations and complement them with 
the matter field equation to complete the evolution system. Numerical considerations 
dictate, however, a preference for the constraint equations (which here act in the outgo- 
ing null hypersurfaces, i.e. they are ODEs in the radial coordinate r). Choosing these 
two hy-persurface equations instead of the other wave-type Einstein equations naturally 
ensures enforcement of the constraints. Unconstrained (or free) evolution systems let 
the solutions drift off the constraint surface of the Einstein equations and may excite nu- 
merical unstable modes, while constrained evolution cuts back on numerical constraint 
violations and offers greater stability. 

There is an analysis of the Einstein equations for the original axissymmetric vacuum 
Bondi problem according to which the equations decompose into the following sets (see 
Refs. [Hbb2] and [dT92b]): 

• four main equations: 

- dynamical equation: 

Roe = (3.2.1) 

- hypersurface equations: 

Rrr = Rr6 = S^'^RaB = 0, (3.2.2) 

where A,B - 6, (p. They contain only derivatives that act in the u — const 
hypersurfaces 

• symmetry conditions: 

Rucp = Rr4, =Re^^O (3.2.3) 
These components vanish identically because of axisymmetry. 

• trivial equation: 

Rur = (3.2.4) 
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• supplementary (conservation) conditions (for energy and angular momentum) 

Rmi = Rue = (3.2.5) 

If they are imposed on a world line, the conservation conditions reduce to regular- 
ity conditions on the vertices of the null cones. They are the analog with respect 
to an r-foliation of the 3+1 momentum constraints. 



It follows by the Bianchi identities that, if the main equations are satisfied, the trivial 
equation is automatically fulfilled and the conservation conditions are satisfied on a 
complete outgoing null-cone if they hold on a single spherical cross-section (e.g. at 
infinity). 

Our situation is, of course, entirely different, we work in spherical symmetry and 
have a matter field coupled to gravity. Still, one might expect to obtain some useful 
insights into the hierarchy of the equations by attempting a similar analysis. We now 
show how Einstein's equations split into hypersurface equations, and equations which 
are automatically satisfied if certain conditions are met. 

In the context of Bondi coordinates, we often abbreviate partial derivatives of a 
function f{u, r) with respect to u and rhy f and /', respectively. 

The hypersurface equations follow indeed from straightforward analogues of equa- 
tions (3.2.2): 

Grr = 8nT,r (3.2.6) 



yields 

while from ^ 



jS' = Inricp'f, (3.2.7) 

/«Rab = 87t/« {tab - IgABT^ , (3.2.8) 
we have, using equation (A.0.21), 

y^Ree = 871 (t - \^^gABT) = 0, (3.2.9) 

and finally ^ 

Y = e^^. (3.2.10) 

The dynamical equation (3.2.1) is satisfied if the hypersurface equations (3.2.7) and 
(3.2.10) hold everywhere and the energy momentum tensor is covariantly divergence 
free: Inserting the hypersurface equations in equation (3.2.1) leads to 

Vr(\,'(\," + r (c/)')^ Y - 2r^(p'(p' + V [cp'f - 2r(p(p' = 0. (3.2.11) 



^We may rewrite Einstein's equations as R„(, = 8n(Ti,i, - ^§"„i,t), where T is the trace of the energy 
momentum tensor of the scalar field. 

^Equation (3.2.10) can also be derived from the (u, r) component of Einstein's equations when plugging 
in equation (3.2.7). 
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This equation is satisfied if V"Tar - 0. 

Clearly, this equation cannot contribute to the dynamics in our setting, since in 
spherical symmetry there are no gravitational degrees of freedom and the dynamics 
resides entirely in the matter field. ^ 

The symmetry conditions, equations (3.2.3), vanish identically because of spherical 
symmetry. What remains to be discussed is the analogue of the trivial equation, equation 
(3.2.4), Gur - SnTur and one of the conservation conditions (3.2.5), Guu = SnTuu- 

The trivial equation is satisfied if the hypersurface equations (3.2.7) and (3.2.10) hold 
everywhere: 

Gur - 8nTur 

^(e^^+2Vp' - V) = 8n^j{(P'f (3.2.12) 

AnrV{(p'f = AnrVicp'f 

Assuming that the hypersurface equations (3.2.7) and (3.2.10) hold everywhere, the 
conservation condition Guu - 8nTuu can be rewritten as 

-2V^ + V ^8nr^i^{(pf -j(p(f)'y (3.2.13) 

At the world line of the central observer, r = 0, this equation is satisfied if the regularity 
conditions, equations (2.3.77), imposed by spherical symmetry hold. The u- component 
of the contracted Bianchi identities VGab - 0, which must also hold for the energy 
momentum tensor as a consequence of the Einstein equations, yields the following 
equation 

V"Euu = -y'Eur, (3.2.14) 

where we use the shorthand E^b = Gab ~ 8nTab- Since we have already shown that 
Eur — everywhere, we have 

= VEuu = -e-^f'VrEuu, (3.2.15) 

and thus 

VrEuu = 0. (3.2.16) 
Combined with Euu = at r = 0, we have shown that Euu = everywhere. 

If one considers the set of hypersurface equations plus the matter field wave equa- 
tion, it is apparent that three integration constants will be necessary to carry out an 
actual time evolution. These constants have been fixed by our choice of proper time at 
the origin. Otherwise, for asymptotically flat coordinates, as in Bondi's original work, 
one chooses the integration constants at infinity. In section 3.6.1 will make an acquain- 
tance with these integration constants (which in our gauge are functions of u) and we 
will see that the conservation condition £;,„ = imposed at future null infinity, J^"*", 
yields a differential relation between two of them, the Bondi mass-loss equation. 

^ Birhoff's theorem shows that a spherically symmetric vacuum solution is necessarily static. The 
conservation of mass and angular momentum in general relativity prohibits the existence of spherically 
symmetric and of dipole waves in the gravitational field. The lowest possible syrmnetry is that of a 
quadrupole which already requires at least axisymmetry. 
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3.3 Field Equations for the Gravitating Massless Scalar Field 



Given the massless scalar (Klein-Gordon) field, all gravitational quantities can be deter- 
mined by integration along the characteristics of the null foliation, since the matter field 
is the only dynamical field in spherical symmetry and needs to be included to make the 
system non-Schwarzschild. This is a coupled problem since the scalar wave equation 
involves the curved space metric. 

The field equations consist of the wave equation for the scalar field 







and two hypersurface equations for the metric functions: 

ti,r = 2nr{(\>,rf 



(3.3.1) 

(3.3.2) 
(3.3.3) 



Note that both jS and V are monotonic in r. Combined with the gauge conditions 
imposed in section 2.3.4 monotonicity implies 



/3 >0 
Vjr > 1. 

In spherical symmetry the curved space d'Alembertian is 



^ + - \dr--du- Idudr + -dry 

r- \r j ,r) r r 



(3.3.4) 
(3.3.5) 



(3.3.6) 



We will now derive a useful relation between the four-dimensional wave operator and 
the wave operator in the two-dimensional (u,r) submanifold. The intrinsic metric in the 
(u,r) submanifold is 

{ds\ = -e^^ {^du^ + Idudr^ (3.3.7) 



and thus the two-dimensional d'Alembertian follows as 



lv\ V 
- dr-2dudr + -dr 



(3.3.8) 



Now we introduce a reseated field which factors out the known falloff of (p at large r 

4> = rep. (3.3.9) 
It is then straightforward to derive the following identity: 



-4, 



(3.3.10) 



*This can easily be obtained by using the formula V^T" = -^d„{ yJ-gT") and setting T" = g^'''Vi,cp. 
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and use this to write the wave equation Dgcf) = as 

Uhxp-y-^ =0. (3.3.11) 

The motivation for introducing the rescaled field i/^ is the following: The amplitude of an 
outgoing spherical wave packet decreases with 1 jr, a rescaled field ip = rep thus behaves 
similarly to a plane wave. Indeed the following identity holds: In the flat background 
case, the field satisfies the usual 1 + 1-dimensional wave equation nfjip = 0. Thanks 
to the plane wave behavior of ip, numerical accuracy is expected to benefit for large r 
since the amplitude of the wave does not change as fast. 

There exist two standard algorithms for solving the wave equation based on this 
identity, which will henceforth be called the Piran-Goldwirth-Garfinkle (PGG) and the 
NSWE (North-South-West-East) "diamond" algorithm which are described in sections 
3.4 and 3.5, respectively. 



3.4 The PGG Evolution Algorithm 

Here we briefly discuss a characteristic evolution scheme due to Piran, Goldwirth and 
Garfinkle [GP87], [Gar95]. We essentially follow the original algorithm although the 
terminology is slightly different. 
It is useful to define the quantity 

f{r,u):^{rcp{r,u)lr (= i/;,,), (3.4.1) 

so that 

f{r,u) - (h{r,u) if, 

(P{r,u), = ^^' \ cp{r,u) = - f{r',u)dr'. (3.4.2) 

^ ^ Jrn 



The wave equation (3.3.11) in w - r coordinates then becomes 



IV 1 (V\ 

/--27^-2(7)/^-^) (3.4.3) 



According to equation (2.3.28), the left hand side can be written as 

df{r{u),r) liV 



(-)^(/-c/)), (3.4.4) 



du 2 \ r 

if the radial positions follow null geodesies. 

The evolution system is thus formed by the ODE's (3.4.1, 3.4.4, 2.3.28). 



3.5 The Gomez- Winicour "Diamond" Algorithm 

The basic idea of this algorithm (see Refs.[GW92b], [GW92c],[Win05]) is to integrate 
the wave equation over the null parallelogram E spanned by the points N, S, W, and 
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E (see figure 3.1 on page 36). Using the identity (3.3.11) and the volume element 
^f-h du Adr = e^^ du A dr We have 

^nh\l> = J dudr(^^^ y (3.5.1) 

E E 

The crucial trick of the whole undertaking is the following: Since h is a two- 
dimensional metric, it is conformally flat. As we will shortly see, □;, has conformal 
weight -2 whereas the surface area element d^x yf—h has conformal weight +2. Thus, 
the surface integral over Djjip is identical to the flat space result which can easily be 
obtained. 

Let gab - ^^gab be a conformal rescaling of the metric gab- First we need to know 
how the determinant of the metric transforms under conformal transformations. For a 
two-dimensional metric ^^^g we have 

det<^'^ cx e^g^,f c ClY'gABf 

and therefore 

i.e. the volume element in two dimensions has conformal weight +2. The formula for 
the d'Alembertian is 

Dg^p-^^a{^P§tM) 
Applying the previous result yields 

V 8 

i.e. the two-dimensional d'Alembertian has conformal weight -2. So, all in all we can 
rewrite the original surface integral in the following manner: 

J d^x ^f-h Dhip = J dudv ^-hfiat □ /m^' 

E E 

The flat two-dimensional (Minkowsi-)metric ds^ = —dt^ + dr^ can be rewritten in 
double-null coordinates as ds^ - -dudv. The metric is {niAB) = ( _i/2 

).So that the 

d'Alembertian becomes Ofjat = -Adudv Now our integral can easily be evaluated: 

dudv V-m dfiatip - -2 dudvip^uv 

E E 

dv{\p^Xi - 4^,v\u0) 

-oO 

= -2(i/'|u1,d1 - ^\ul,vO - ^\uO,vl + 4'\u0,y0) 

= -2i^iN)-ipiW)-ipiE) + ipiS)) 



^See [Wal84] Appendix D for a definition. 
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where we have introduced the grid points N, S, E, and W as shown in figure 3.1 



Figure 3.1: The computational molecule for the NSWE scheme is made up by two u = const 
hjrpersurfaces and two ingoing v = const null geodesies. 



3.6 Compactification of the Radial Coordinate 

Compactification is rather simple for characteristic codes, and has been used extensively 
in the characteristic approach to numerical relativity [GW92b, GW92a, GW92c, WinOS]. 

The type of compactification we pursue here is not conformal compactification per 
se (see Ref. [riau4] for a review), where the spacetime manifold is completed by 
"attaching" a boundary at "infinity". The method of conformal compactification maps 
physical spacetime onto a bounded open region of unphysical spacetime, introduces an 
unphysical metric via a conformal rescaling g^b - Cl~^gab (O, approaches zero (at infinity) 
at an appropriate rate), and factors out known asymptotic behavior of geometrical 
quantities by a conformal transformation. 

Rather, we use an ad-hoc regularization adapted to the simple geometry which 
could be related to conformal rescalings [Hus07]. We employ a simple coordinate 
transformation with respect to the radial coordinate, which maps a half-infinite domain 
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Using (3.3.11) and the volume element 
written as 





(3.5.2) 



r e [0, oo) to a compact interval x £ [0, 1]. In this ad-hoc approach the metric is not altered 
(no unphysical metric is introduced) and the evolution equations "notice where infinity 
is" because they degenerate there [Fra04]. The advantage of this compactification 
method are its simplicity and that we can evolve till future null infinity, J^^ and thus 
study global and radiative properties of collapse spacetimes in detail. 

We would like to emphasize that the use of conformal rescalings does not imply 
that one has to employ the rather difficult framework of Friedrich's regular conformal 
Einstein equations [FriSl]. In fact such an approach has also been suggested in [And02, 
HSVZ05]. 

Before introducing compactification we want to say a few words about asymptotic 
series expansions. 

3.6.1 Asymptotic Series Expansions for the Massless Scalar Field 

Assuming initial data that are smooth at J^^ , one can expand the massless scalar field 
(/) in powers of 1/r near J^"*" ^ 

c/)(u,r) = ^ + ^+0(r-^ (3.6.1) 

The coefficient Cnp of the term in the expansion is a Newman-Penrose constant 
[NP68] of the scalar field. 

Inserting the expansion (3.6.1) into the hypersurface equations (3.3.2) and (3.3.3) 
yields 

and 

V{u, r) = e^C) (r - 2M(m) + J + 0{r-^), (3.6.3) 

where integration constants H{u) and M{u) have been introduced. 

H{u) - limr^oo fiiu, r)\u=const indicates redshift since Bondi time Ug is related to proper 
time at the center via the relation (2.3.35). 

M{u) = lim |^r2e-2^^(") (Y-^ J . (3.6.4) 

is the Bondi mass which is in general not conserved. In terms of the Misner-Sharp 
mass-function defined in section 3.7.1 we have the relation 

M(m) = \imm{u,r)\i,^const, (3.6.5) 

r — >oo 

and the asymptotic expansion 

m{u, r) = M{u) - + Oir^). (3.6.6) 



^{u, r) = H{u) ^ + 0{r-% (3.6.2) 



*'A constant term in the expansion of (p has been omitted since one can trivially rescale the coupled 
Einstein massless scalar field system: A transformation (p —> cp + const leaves nip = invariant. 
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3.6.2 Compactified Evolution Scheme 



First, we introduce a compactified radial coordinate, then we recalculate our evolution 
equations and aim to regularize them at null infinity. 

Starting with our standard Bondi-like coordinates {u,r,0,(p) we introduce a compact- 
ified radial coordinate 

-=T^ (3.6.7) 

which maps r £ [0, oo] i-> x e [0, 1], so that points at J^'^ are automatically included in 
the grid at x = 1. Consequently, we have the relations 

r = dr = {l- x)-^dx, and = (1 - x)^-^ . (3.6.8) 

1 - X \ J ' dr ^ ' dx ^ ' 

The line element, equation (2.3.17), then becomes 

ds' = -,2^du(y(i^)du + (YT^^^) + (t^J^^'' (3.6.9) 

which evidently contains singular terms. Regularity, however, only requires that the 
field equation for the scalar field and the hypersurface equations for the metric functions 
be well behaved. 

A naive approach of rewriting the hypersurface equations in terms of the x-coordinate 
leads to a singular equation for the quantity V 

ti,, = 2nxil-x){(P,,f, (3.6.10) 
= n ^' (3.6.11) 

(1 - X)2 

while the equation for ji is fine. 

The introduction of a renormalized quantity V :- y(l - x) (chosen such that 7 = 7) 
lets us rewrite the hypersurface equation (3.6.11) in the following form 

V\ e^^ - VIx 

' ' (3.6.12) 



xj^ x(l - x) 

which remedies the situation somewhat. Still, it is not obvious that (7),!- is well be- 
haved at , since the denominator of equation (3.6.12) tends to zero. An asymptotic 
expansion shows that {^),x is nonsingular, 

=2M(w)e2^("U<9(r-i), (3.6.13) 

but the cancellation in the numerator of (3.6.12) that brings about this result is very 
delicate in numerical terms. 
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As described in [t^H AOS], to obtain a fully regular system of evolution equations, 
we eliminate V (or V) by the Misner-Sharp mass-function 



m{u, X) 



which satisfies 



2(1 X 



The set of hypersurface equations then becomes: 



(3.6.14) 



(3.6.15) 



= lnx{\ - x)((/),,)2. 



1 m 



(3.6.16) 



These equations are now completely regular. Note that the term m)-^ = mfr does not 
cause problems because of the smoothness of the metric at the regular center (2.3.77). 

In section 4.1.1, we will choose our gridpoints to freely fall along ingoing radial null 
geodesies x{u) which fulfill 



dx 
du 



= -l(l-.)V/^(l-2^1^). 



(3.6.17) 



In section 4.1.3 we will argue that this choice is crucial to resolve DSS phenomena. 
For ip the matter field equation takes the form 



18) 



At this reduces to the ODE 



(3.6.19) 



The diamond scheme can be derived by applying the transformation to the x- 
coordinate to equation (3.5.2), which yields 



4>N = 4'W + 4'E 



-^s-'^ J dudx(^^-^^2mi 



(3.6.20) 



Equations (3.6.16) and (3.6.20) now form a manifestly nonsingular set of evolution 
equations for the massless scalar field coupled to gravity. 
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The gauge and regularity conditions (2.3.77) outlined in section 2.3.4 and the reg- 
ularity of the scalar field (p (2.3.83) along with the definition of the rescaled field 
equation (3.3.9) now become 

ti{u,x)^0{x% 

m{u,x)=0{x^), (3.6.21) 
ip{ii,x) - 0{x). 

In section 4.1.2 we discuss Taylor expansions in the vicinity of the origin. 



3.7 Diagnostics 

3.7,1 The Misner-Sharp Mass 

In spherical symmetry, there exists a well defined notion of quasilocal energy, the 
Misner-Sharp mass-function (see Refs. [MS64] and [Hay96]): 

m-'^[l-n, (3.7.1) 



which in outgoing Bondi coordinates yields 

(3.7.2) 



m{u,r) = ^ 



Note that m/r is a smooth function. The Misner-Sharp mass measures the energy content 
of a sphere of radius r and reduces the ADM and Bondi masses in the appropriate limits. 
We also refer to the Misner-Sharp mass-function defined in equation (3.7.2) as niMs, in 
order to distinguish it from an integral expression for the local mass, nip which we 
define below. 

We write equation (3.7.2) as an integral along outgoing radial null geodesies {u = 
const, 9 - const, (p - const): 



mp{u,r) = I dfm'{u,f) (3.7.3) 





and then use Einstein's equations to simplify the integrand: 

1 I. V _,a\ r V _,V\ 



(3.7.4) 



Since m' is a form of local energy density, we refer to its integral nip (defined by equa- 
tion (3.7.3)) as the "integrated-scalar-field mass-function". 

On the continuum level, the mass-functions wJms and nip are clearly identical, but 
numerically they will in general differ by a small amount due to having different finite 
differencing errors; their difference is a good "diagnostic" of the code's accuracy. 
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To this end, we define 



bm{u,r) = (3.7.5) 

"^total, init 

where ?«totai,init = '«Ms(w=0/''max) is the total mass of our initial slice, i.e. the mass- 
function at the outer grid boundary on the initial slice, bin is a dimensionless measure 
of how well our field variables satisfy the Einstein equations; we must have \bm\ «; 1 
everywhere in the grid at all times for our results to be trustworthy. 



3.7.2 The Bondi Mass and the News Function 

Taking the limit of the Misner-Sharp mass-function, m{u, r), as r — > co at constant u, we 
obtain the Bondi mass: 

M(m) = limm(w,r). (3.7.6) 

r — >oo 

In an isolated system outgoing waves can radiate physical energy to future null 
infinity, .y* . Therefore, the Bondi mass is in general not conserved in retarded time. 
(In contrast, the ADM mass, which is defined by taking the limit of m{u, r) as r — > oo at 
constant t, that is, on a spatial slice, is conserved in t.) 

Moreover, one can show (see Ref. [Wal84]) that there exists a flux / such that 



AM = - 




(3.7.7) 



where ^ is a cross-section of --^^ with a u - const null hypersurface and 7 is a real 
interval. 

We will derive a relation between the outgoing radiation flux, which is described 
in terms of the scalar news function, and the change of the Bondi mass in time. This 
relation is known as the Bondi mass-loss equation. (As always, u denotes central time and 
not Bondi time. The two time coordinates are related by equation (2.3.35).) 

We insert the asymptotic series expansions of the fields (p, ^ and V, given in equations 
(3.6.1), (3.6.2) and (3.6.3), respectively, into the (w, M)-component of the Einstein equations 

= r2£„„ = (G,„ - SnTuu) 

= + y [^^ -V - 2r^) + V-8n (r(pf + AnVr(^j{cj)')^ - 2^cj)'^ (3.7.8) 

and obtain 

= (0(r-2)) 

+ l^e^ _ - - Irti) + IHe^'' (r - 2M) + O(r-i)) 

+ (e2«(-2M + 6)(r-i))) (3-7.9) 

-{8n{cf+0{r-^)) 

+ (0(r2)(0(r-4)-6)(r-3))) 
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As it turns out, the third and the fourth terms are crucial, all the other terms are at least 
0(r"^) or cancel out. We find 

= lim (r^Euu) = -2e^^^"^M{u) - 8n {c{u)f (3.7.10) 
Now, we introduce the scalar news function as 

N{uc) = e-^^^"c^c{uc). (3.7.11) 
Combining the last two equations yields the Bondi mass-loss equation 

^-2Hiuc)dMiUc) ^ _4„^2(„^)_ (3 712) 

duc 

In Bondi time, it becomes, using relation (2.3.35), 

^ = -AuNHu.). (3.7.13) 

From equation (3.7.7) we find 



JAllQ 

= -471 I N^{u^)du^ ^ - I dCl I dii^N\u^) 



(3.7.14) 



Thus, the square of the news function is just the flux that appears in the integral (3.7.7): 

N\u^) = f. (3.7.15) 

The positivity of the flux-function / entails that 

AM < (3.7.16) 

always. "News", that is radiation, can only decrease the Bondi mass contained on a 
collection of null-slices that extend to But the Bondi mass can never increase to the 
future. 



3.7.3 The Bondi Mass as a Linkage Integral 

Let be a generator of an asymptotic time-translation and let S^a be a one-parameter 
family of spheres which approach the cross-section of J^'^. Then the Bondi mass is 
defined as (see Refs. [GW81] and [Wal84]) 

MBondi = - l^rn -^ f Cabcd^'^'', (3.7.17) 
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where the gauge condition Va^'^ = must be fulfilled in a neighborhood of . For the 
Bondi metric, equation (2.3.17), the condition S/a^'^ - gives 

+ Cr + 2^,ue + 2^,E: = (3.7.18) 

if = ((^",(^',0,0). The generator of (asymptotic) time translations 

^« = (e-2/!^ 0,0,0) (3.7.19) 

satisfies the vanishing divergence condition as can be verified by making use of the 
asymptotic expansions for |S and V given in section 3.6.1. Now the integrand of equation 
(3.7.17) is 

Cam^'E,'^ dx"dx^ = V=^(V".^'' - Ve) de d(p = -2M{u) sin 6 dO d(p (3.7.20) 
so that one obtains 

Msondi = M(m), (3.7.21) 
where M(m) is the 1/r coefficient in the asymptotic expansion of ^^"''"^ . 

3.7.4 The Ricci Scalar Curvature 

From the Einstein equations (3.1.3) for a minimally coupled massless scalar field (p one 
finds by taking the trace 

tr Gah = -R = -871 tr T^b = -87Z(/)'"(/),„ (3.7.22) 
the following expression for the Ricci scalar 

R = 87z/''(/),«(/),fc. (3.7.23) 
For the Bondi metric (2.3.17) this yields 

R = 87ze-2^*(/)'(-2(/) + jcj)'). (3.7.24) 

3.8 Double Null Equations 

To complement the DICE code which is based on Bondi coordinates, we introduce a first 
order formulation of the Einstein-scalar field system in double null coordinates which 
offers the advantage that it can penetrate apparent horizons. 

We consider the Einstein massless scalar field equations Rah - 8n(pa4>b for the spher- 
ically symmetric metric introduced in section 2.4 

ds^ = -a\u, v)dudv + r^{u, v)dd?- (3.8.1) 
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in double-null coordinates {u,v). The tensor components are given in appendix B. We 
often use shorthands for partial derivatives: 



d_ 
du 



dv' 



(3.8.2) 



We set s := '\4n(p to simplify the appearance of the matter field terms and define 
additional variables in order to write the equations in first order form following Ref. 
[HS96]. 

First, we define the following evolution variables 

(3.8.3) 
(3.8.4) 
(3.8.5) 
(3.8.6) 



Dl: 


P 


:= s 


D2: 




■=s' 


D3: 


f 


._ y. 


D4: 


g 


■=r' 


D5: 


d 


._ «' 

a 



(3.8.7) 



The Einstein and matter field equations then become 



El: 


/' 


E2: 


d 


CI: 


f 


C2: 


g' 


SI: 


4 


S2: 


P' 



-pq 



r 

f(] + gP 



(3.8.8) 

(3.8.9) 

(3.8.10) 
(3.8.11) 
(3.8.12) 

(3.8.13) 



where we denote Einstein equations of wave-type by the letter E, constraint equations 
by the letter C, and the scalar field wave equation by the letter S. 

Boundary conditions at the center of spherical symmetry are dictated by regularity 
and gauge choices as detailed in sections 2.4.2, 3.8.1 and 3.8.2 and 



ds 
dr 
da 



g = -f = a/2 

= 
= 



(3.8.14) 
(3.8.15) 

(3.8.16) 
(3.8.17) 
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To evolve the coupled first order system, we choose a constrained evolution using 
equations E2, SI and D2, D5, C2, D4. This leaves equations CI, Dl and D3 for checking 
the numerical solution. 

We need to specify the variables a, r, s, p, f, g at the axis, while q and d are calculated 
from evolution equations. Obviously, r = at the axis. 



3.8.1 Regularity and Boundary Conditions 
Regular Variables at the Center 

We demand that our evolution variables (metric functions, scalar field) be well defined 

at the center of spherical symmetry. See [ ] for a detailed analysis for Cauchy 

problems. As we have seen in section 2.3.4, this is most easily obtained by demanding 
that the variables must not have a kink as the origin is crossed along a ray parametrized 
by linear distance /. Effectively, this forces all odd derivatives of the metric with respect 
to r at constant f to vanish at the origin. 

Therefore, a^{u, v) must be an even function of r, i.e. 



so that we arrive at the boundary condition 

da 
dr 



0. 



ds 
dr 



Similarly, we have for the scalar field 

0. 

We can also investigate the behavior of g,-,- - - 

grr = g%+0{r') 



(3.8.18) 



(3.8.19) 



(3.8.20) 



(3.8.21) 



Local Flatness at the Center 

As it turns out, one has to also impose local flatness to ensure a regular behavior in the 
evolution equations. Again, see [ aG05] for details. Local flatness means that the spatial 
metric locally looks like the flat metric 

df = dr^ + r^dOt^. (3.8.22) 

In double-null coordinates (w, v) the flat spacetime metric is given by 

ds^ = -a^dudv + r^dQ^ (3.8.23) 
^We may use flat space coordinates (t, r) since we impose local flatness later. 
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where a - a{u,u) - const and r is given by 



r="-^. (3.8.24) 
It then follows that r' — -r or, in the notation of section 3.8), 

g--f-\ (3.8.25) 

and also that fg + a^/4: = 0. 

Due to spherical symmetry and local flatness, r has to decrease as fast (with u) on an 
ingoing (radial) null ray as it has to increase (with v) on an outgoing null ray close to the 
origin. Also, a null ray through the origin only appears to change direction in spherical 
coordinates, while it obviously passes straight through in Cartesian coordinates. 

As we impose local flatness at the center, we may use flat space {t, r) coordinates to 
express boundary conditions, such as ^,\t - 0. It is straight forward to implement such 
boundary conditions on an equispaced double null grid. 

Issues of regularity of the evolution equations are discussed in section 3.8.2. 

3.8.2 Regularity of Evolution Equations 

At r = 0, the formally singular right hand sides of equations El, E2, and SI, S2 yield the 
following conditions: 

fg + fl2/4 ^ (3.8.26) 
dr(fg + a^/A) = (3.8.27) 
fq + gp = (3.8.28) 

The last condition leads to the final boundary condition needed, which is simply 

p^q. (3.8.29) 

3.8.3 Diagnostics for the Double-Null System 

We calculate the density function 2m/r, 

Imir = 1 + (3.8.30) 
where m is the Misner-Sharp mass, defined in section 3.7.1, 



(3.8.31) 



The quantity maxi, 2m jr is used as an indicator for critical behavior and for the closeness 
of a slice to the formation of an apparent horizon. To detect and measure an apparent 
horizon we search for a zero in the function g = r-o (see 2.4.1). 
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The scalar curvature can be expressed in terms of the scalar field or the geometry: 



R = -8pq/a^ 

2a^ + 8r^ {ad' - da') + Sa^ {2rf' + fr') 



(3.8.32) 



The scalar field energy density is given by 



P = 



(3.8.33) 



It is instructive to introduce observers at constant radius r{u,v). We would then 
like to compute the proper time t{u,r) along each of these worldlines of constant r. 
Following [HS96] we invert r = r{u,v) to obtain v = v{u,r). At r = const, it follows that 

dv I 



dv - ^\rdu. The line element (2.4.1) then becomes 



ds \ = -a (u,v) — 

\r=const ^ ' 



dip- + r'^{u,v)dQ^ 



Then, we have 



r-tl r-u 

t{u,r = const) = ( ^|-Suudu = I a{u,v{u,r)) ^ — dii 
Jo Jo V C^M r 



/o JO 
At r = 0, local flatness implies v{ii, r) - + u, and thus 



t{u,0) 



r 

I a{u,u)du 
Jo 



(3.8.34) 



(3.8.35) 



(3.8.36) 
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Chapter T 

Numerical Algorithms 



4.1 The DICE-code 

The compactified code used in [PHA05] is based on the "DICE" (Diamond Integral Char- 
acteristic Evolution) code, which has been documented in [HLP"'"00] (there particular 
emphasis is given to detailed convergence tests). 

4.1.1 Our Overall Computational Scheme 

While this section discusses the uncompactified Bondi evolution algorithm, it can be 
applied to the compactified system, by a few simple replacements. Most notably, the 
metric function V has been replaced by the Misner-Sharp mass m. See section 3.6 for 
further details. 

Summary 

We first construct initial data: 

1. Choose the field i/^(mo/ ^) on a null slice u = uq. 

2. Also choose the positions, rj(uo), of our grid points on this same initial slice. 

3. Radially integrate the hypersurface geometry equations (3.3.2) and (3.3.3) out- 
wards from the origin to compute jS and V/r at all the grid points on the initial 
slice. 

Our integration scheme starts by first integrating the geodesies back in time one step 
Am, since our geodesic integrator needs two slices to work. Now suppose we know all 
the gravitation and matter fields for time levels < k. To determine them for time level 
k + 1, we use the following algorithm (we use the usual notation where superscripts 
denote time levels): 
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1. Determine the Taylor expansion of ip (r) near the origin (cf. section 4.1.2). 

2. For each grid point in some small (typically 3 grid points) "Taylor series region" 
starting at the origin and working outwards, 

(a) Integrate the geodesic equation (4.1.3) ahead one time level to determine the 
position of this grid point at time level k + 1. 

(b) Using the ip^ Taylor coefficients and the Einstein equations (cf. section 4.1.2), 
compute first ijj'^'^^, then (i^*^ and (V/r)'^*^ at this grid point 

3. Now for each grid point in the rest of the grid, starting just outside the Taylor 
series region and working outwards, 

(a) Integrate the ingoing null-geodesic equation (2.3.28) ahead one time level to 
determine the position of this grid point at time level k + 1. 

(b) Use the Winicour diamond scheme (3.5.2) to compute ip^'^^ at this grid point. 

(c) Radially integrate the hypersurface geometry equations (3.3.2) and (3.3.3) 
outwards one spatial grid point to compute f5^*^ and (V/r)'^*^ at this grid 
point. 



Discretization of the Evolution System 

We evaluate the integral in equation (3.5.2) by treating the integrand as constant over 
the (small) null parallelogram L, as shown in figure 3.1. Thus we can compute the 
integrand at the center of L which in turn equals its average between the points W and 
E to second order accuracy: 

Integrand at center = h 0{h ), (4.1.1) 

where h = ^ yjiAu)^ + (Ar)^. Then we multiply the integrand by the area of the null 
parallelogram E, which can be approximated by 

As will be discussed in section C.3, this algorithm is globally 2nd order accurate; for 
the special case of a flat background ((7)^ = 0) the algorithm is exact. 

We use the explicit trapezoidal method (see equation (C.2.8) in appendix C) to 
discretize the hypersurface and ingoing null geodesic equations. In principle the dis- 
cretization is straightforward, with the exception of the null geodesic equation. The 
corrector requires evaluating the right hand side at the n + 1 time-level, but this has not 
been computed yet at the time we do the geodesic integration. To solve this problem, 
we radially extrapolate the needed V/r value from V/r and (V/r)' one spatial gridpoint 
inwards at the same n + 1 time-level. At the origin, we know V/r — 1 thanks to local 
flatness, an}rway. 
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Freely Falling Gridpoints 



The original Winicour algorithm uses a fixed grid (i.e. the radial positions of the grid 
points do not move in the course of evolution. As the coordinate velocity of light is 
in general not constant in a curved space-time, the vertices N, E, S and W of the null 
parallelogram E cannot be chosen to lie exactly on the grid. Therefore one has to use 
interpolation to compute the matter field and metric functions at the vertices. 

In the DICE code we used a different approach (which is, in fact, borrowed from a 
different characteristic evolution scheme due to Piran, Goldwirth and Garfinkle [ ] 
as described in section 3.4): Our gridpoints free-fall along ingoing null geodesies. To 
compute the locations of the gridpoints on a new slice we integrate the ingoing null 
geodesic equation 

dr_ 
du 



1 V 

-- . (4.1.3) 



With this approach, the vertices of the null parallelogram E in the diamond algorithm 
(cf. section 3.5) are gridpoints, so we don't have to interpolate there. 

As well as avoiding interpolation, this scheme has the additional advantage that it 
provides a sort of adaptive grid refinement for free, since, in regions of strong gravity, 
the null geodesies are being focussed and thus the density of gridpoints is increased. 

Since our grid points move along ingoing null geodesies, each grid point eventually 
reaches the origin. At this point it disappears, i.e. we remove it from the grid. 



Stability of the Winicour algorithm and adaptive choice of timesteps 

Assuming that the gridpoints follow ingoing null geodesies, we want to impose a 
restriction on the "drift" of individual gridpoints from one slice to the next. It follows 
from the null geodesic equation (4.1.3) that, if we want the change in r from a slice u" to 
u""*"^ to be smaller than a factor D ('driftJimit') times Ar, we have to demand 

All <D'^^r (4.1.4) 

for all gridpoints. (The minus sign in the null geodesic equation is absorbed by the Ar 
which is a directed quantity.) The code assumes that D <= 1, i.e. we drop at most one 
grid point per time step. 



4.1.2 Taylor Expansions 

To start off the integration of the evolution system on a new null slice at retarded time 
u = iiQ + All and time-level k + 1, we Taylor expand the scalar field on previous slices, 
and employ the matter field wave equation to propagate the information to the new 
slice. 
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Calculation of ip 



First, we compute a Taylor expar\sior\ of the rescaled matter field ip on the slice at u = uq 
and time-level A: in a region that encompasses the origin of spherical symmetry, r — 0: 



ipHuo,r) = tor + hr^+0{r^). 



(4.1.5) 



We have omitted the constant term due to the regularity of the matter field (j) and the 
definition of the rescaled field ip — cpr. For the sake of simplicity, we determine the 
coefficients Iq and ti by fitting a parabola through two sample points (and the origin). 
A fancier approach used in the code is to use a general linear least squares fit over a 
number of gridpoints adjacent to the origin, say 5 points. This way of handling the 
origin has in practice proven valuable in smoothing out discontinuities in the fields 
near the center. 

Next, we establish a relation between the field values on the new slice, at time-level 
k + \, and the field on the old slice, at time-level k. To achieve this we make use of the 
field equation 

'V\ e-^Pip 



= 0, 



where 



Solving for the ^,,-term we arrive at 



(7)' 



'V\ V 

dr - Idudr -\ drr 

r 



V\' /il> \ V 



(4.1.6) 



(4.1.7) 



(4.1.8) 



Replacing the u derivative by a finite difference approximation 



{duff 



Am 



+ (9(Am) 



yields 



Aw 



Using the behavior of V near the origin 

V , 
- = 1 + 0{r^) 



r 



and 



(4.1.9) 



(4.1.10) 



(4.1.11) 



and integrating the above relation, equation (4.1.10), with respect to r we recover 

Aw 



i/;'^+i(r) = i/;«(r) + 



dr[r^(T)+0(T)\. 



(4.1.12) 
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Finally, we insert the Taylor expansion for ip into equation (4.1.12) and find the desired 
result 

xp'^+^r) = tor + hr^ + h^u r + 0{^^), (4.1.13) 

where A denotes a small quantity (either r or Aw). 

The calculation works the same way for higher orders. To third-order accuracy we 
find 

xp^*^{r) = Iqt + hr^ + firAw + '^t2r^M + ^t2r{Auf + t2r^ + (9(A^). (4.1.14) 
Calculation of the Metric Functions 

We can determine the metric functions on the new slice by inserting the expansion of 
the rescaled field, equation (4.1.14), into the hypersurface equations 

P,r^2nr(cl),rf (4.1.15) 
y , = e^l^. (4.1.16) 

(P(r) ^to + tir + tiAu + 0{A^) (4.1.17) 

^'{r) = 2nr{h)^ +0{A^) (4.1.18) 
p{r) = n {hf + 0{A\ (4.1.19) 



We first determine 



and thus 
and 

Along the same lines we have 



Y{r) = 1 + 2|g + 0{<^^) = 1 + 271 {hf ? + 0{^) (4.1.20) 

and eventually 

\r)^\^\n{hfr'^^0{A^) (4.1.21) 

(y)'(r) = ^7z(ii)2r + 0(A2). (4.1.22) 
The third-order accurate expansions are 

jg'^+i(r) = 87Z Qf2,2 + + \hhr'] + 0{A% (4.1.23) 

and 

y'^+i(r) = r - + 871 (^^^^) + C)(A^). (4.1.24) 

The Ricci Scalar curvature (see section 3.7.4) has the following expansion at the 
center: 

R{u,r = 0) = IGntj. (4.1.25) 
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Taylor Series for the Compactified Code 



For small r the functional behavior of the x-coordinate is similar to r, so that we can 
reuse the Taylor series expansions of the fields as functions of (u, r) as given in section 
4.1.2 and transform them to the x-coordinate by applying the expansion 

r= ^^^=x + x^+x^ +0{x''). (4.1.26) 

To third order accuracy we obtain (A denotes a small quantity Am or x) 

+ {to + 2h + t2) x^ + 0{A'^), (4.1.27) 

^^*^{u,x) = (nt\ + 3nht2Au)x^ + {lnt\ + ^^1^2)^^ + 0{A^), (4.1.28) 
We include one more order for m, since it appears as m/x in the field equations. 



m 



\u,x) = i^-^t\ + 2ntit2Au - -^t\Au^^x^ 

(1 9tl \ 

- + lnt\ + lnht2 - nil Am - ^"^ ) + 0{A^), (4.1.29) 



and 



(^) ("'^^ ^ {lnt\ + 6ntit2Au - '^t\Av?-^x^ 

+ (2 + %nt\ + Sntit2 - ^ntj Am - 9ntl Au^) x^ + 0{A'^). (4.1.30) 



4.1.3 Mesh Refinement 

Hamade and Stewart [HS96] have implemented full Berger-Oliger mesh refinement in 
double-null coordinates (without compactification) to achieve sufficient resolution to 
study critical collapse. Garfinkle [Gar95] has shown, that this is not really necessary - 
here we follow his approach to increase resolution: Most importantly we choose our 
gridpoints to follow ingoing radial nullgeodesics. This leads to a rapid loss of gridpoints 
in the early phase of collapse, but to an accumulation of gridpoints in the region of strong 
curvature for the late stages of critical collapse (see figure 4.1). Furthermore, when half 
of the gridpoints have reached the origin we refine the grid and thus obtain a very simple 
but effective form of mesh refinement which is a crucial ingredient in the calculation 
of critical collapse spacetimes. In previous work [HLP+00] we have also tuned our 
outermost gridpoint to be located just outside of the past self -similarity horizon (SSH), 
the backwards lightcone of the accumulation point of the self-similar solution. Here 
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Figure 4.2: The convergence of the error diagnostic Euur with increasing grid resolution for 
two near-critical evolutions. Evolution (1) uses 5000 gridpoints and p = p*[5000] + 10"^" and 
evolution (2) uses 10000 gridpoints and p = p*[10000] + 10"^''. The data displayed as dots and 
circles have been sampled. 



we choose to go out all the way to null infinity. The most effective approach in this 
situation would be to just refine the region inside the SSH when half of the gridpoints 
in this region have reached the origin. While this is straightforward to implement, we 
found the penalty on the resolution that the original condition (the loss of half of all 
gridpoints) causes to be acceptable for the results presented here. 



4.1.4 Accuracy and Convergence 

The code is globally second order accurate. Figure 4.2 shows a convergence test for near- 
critical evolutions. We want to emphasize that the critical value p* of the initial data 
parameter depends on the grid resolution. This fact is essential when doing convergence 
tests for near critical evolutions, as has been discussed in our previous paper [HI.,P+00]. 

To monitor the accuracy of the code during runs we use components of the Einstein 
equations which are automatically satisfied if the evolution equations hold, such as the 
following linear combination of the (u, u) and (u, r) components 

- E„„, = i^(Guu - SnTuu) - rHv/r)(Gur - SnTur). (4.1.31) 
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Figure 4.3: A 3-level convergence test for ip{r) 



This can be rewritten in the following form 

Euur = 2e^l^m + 8n 

where / = Since this expression is a linear combination of tensor components, we 
use a suitably normalized quantity Euu,- = TTEj^' where E|Mt,r| is the sum of the running 
maxima (max,<y over au = consi-slice) of the absolute values of the individual terms of 
Euur- We must have Euur 1 for our finite difference solution to be a good approximation 
to the continuum solution. A convergence test for this quantity is shown in figure 4.2. 

We also display numerical convergence tests (see appendix D) for ip and m (3- 
level) and 6m (2-level) in figures 4.3, 4.4, and 4.5 which show that our code computes 
these gridfunctions with the accuracy expected from the theoretical considerations of 
appendix C.3. 



4'- 



12 



(4.1.32) 
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Figure 4.4: A 3-level convergence test for niysir)- 
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Figure 4.5: A 2-level convergence test for 6m. 
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4.2 The Double-Null Code 



The numerical solution scheme for the double-null code follows the work of Hamade 
and Stewart [ ] and further improvements made by Harada and Carr [Har05] which 
makes the scheme second order accurate. 

We choose to refer to the evolution equations by the mnemonic abbreviations Dx, 
Ex, Cx and Sx, (with x, a number), given in section 3.8, rather than usual equation 
numbers. 



4.2.1 Regularization 

At the center of spherical symmetry we need to regularize equations El, E2, SI and S2 
using I'Hospital's rule. To simplify matters we use derivatives evaluated at r = 0. 
Derivatives of the type are straightforward to write down on an equidistant grid in 
u and V. See figure 4.6 for an illustration. However, we need finite difference approxi- 
mations for unevenly spaced grids, as our grid is only evenly spaced in u and v, but not 
in r. Such finite difference approximations can be obtained by determining a low order 
polynomial through the required number of points and using symbolic differentiation. 
This is most easily done in a computer algebra system, such as Mathematica [WR05]: 

For example, a second order accurate approximation of the first derivative of a 
function /, given the function values (/o,/i//2) at three points {roJiJi), evaluated at 
the point yq, given by 

/(r2)(r0 - rl)^ + (2r0 - rl - r2)(rl - r2)/(r0) - (rO - r2)V(rl) 
^ ^''^ = (r0-rl)(r0-r2)(rl-r2) ' (^•'•'^ 

can be obtained by the compact expression 

rgrid = {r0, rl, r2}; 

FullSimplify [NDSolve ' FiniteDif f erenceDerivative [ \ 
Derivative [1] , rgrid, Map[f, rgrid], \ 
DifferenceOrder -> 2]][[1]] 

Similar expressions can be derived for higher order accuracy and second derivatives. 

Depending on the number of previous slices that we have already computed, we 
employ different general finite difference approximations for these radial derivatives at 
constant f. As it is not always possible to use second order accurate approximations, 
we have to resort to lower orders in the first two timesteps of an evolution. This is not 
a big problem, though, if we stick to initial data of compact support, such as Gaussian 
data, where the field initially almost vanishes at the origin. 

To regularize the E2 equation d = -^-^ pq according to I'Hospital's rule, we 

take the second derivative of the numerator of the singular part of the right hand side, 
fg+a^/A, using a second order accurate general finite difference stencil that requires four 
points. All other singular equations require only the computation of first derivatives as 
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r = 



Figure 4.6: This figure shows the double-null grid around the origin. We have indicated the 
new slice being computed by (with time index +) and the v-index that hits the origin on that 
slice, by /O"^. Four gridpoints which lie at constant t (flat space time coordinate) are indicated by 
filled circles. Using these gridpoints, one can approximate derivatives of type ^ |f . 
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the denominator of the singular right hand sides are r. To regularize equation SI we 
compute an approximation to - {^fg + using formula (4.2.1), given above. 

The right hand side of equation S2 is identical to that of equation SI, so that we can 
reuse the right hand side approximation computed for SI. As it turns out, the right 
hand side of equation El is in fact zero analytically. This is due to regularity of equation 
E2 which, in addition to fg + = 0, also forces 



difg + a^/2) 



dr 



(4.2.2) 



4.2.2 Initial Data 



We may freely specify unconstrained data for the scalar field s{u = const, v) on an 
initial nullcone, usually chosen to lie at u = 0. For the sake of simplicity, we restrict 
considerations to Gaussian initial data: 



s(0, v) = A exp 



(4.2.3) 



Together with the specification of the gauge for the null coordinate v, such as 

ri0,v) = v/2, (4.2.4) 

we can construct initial data for all remaining evolution variables, a,d,f,g,p, and q. ^ It 
is then trivial to reconstruct c]{0,v) and g{0,v) by analytic differentiation and also easy 
to compute d{0, v): 

D2: qiO,v) = s'(0,t;) = -2^^^s(0,z;) (4.2.5) 

D4: ^(0,z;) = / = 1/2 (4.2.6) 

2 

C2: dM = ^. (4.2.7) 

The remaining quantities p{0,v), f{0,v) and a{0,v) have to be reconstructed by nu- 
merical integration: 

D5: (lnfl(0,i;))' = - =d (4.2.8) 

a 

fg + a^/A 

El: f'{0,v) = - ^^ ^ (4.2.9) 

S2: p'{0,v)^-^^-^. (4.2.10) 

We use Simpson's rule (C.2.12) (and the Trapezoidal rule (C.2.11) for the first gridpoint 
outwards from the origin) to perform the integration. 



^Note that Hamade and Stewart [ ] fix the gauge by choosing d(0, v) = 0, instead. 
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We also apply I'Hospital's rule for formally singular right-hand sides at the center 
as described in section 4.2.1. Only the equation for p is formally singular at the origin. 
Since we do not know any previous slices, we cannot use derivatives and instead 
use -^lu and the gauge choice (4.2.4): 



p'\o = - 



-2 



dv 



u,0 



u,0 



if(l + gP) 
{fq+p/2) 



(4.2.11) 



= -2(fqy \o-gp'\o 

Subscripts of a gridfunction denote v indices on the numerical grid, starting at / = 
at the origin. Using this expression we can derive the following approximation for pi 
by applying the Trapezoidal rule 

pi = 



Av 



Ml 



) 



1 4- ("^ 4- Sl\ 

^ ^ 2 \Av ^ ri) 

Similarly, we can derive an expression for p2 by applying the Simpson rule 

PO + 2Av ^ IK^ 4 



(4.2.12) 



Pi 



For indices i > 2 we have 



1 4- 42 (-S0_ , g2\ 

3 \ Av^ Vi) 



(4.2.13) 



P,-2 + -3- 



Pi = 



^ fl+gp 








i-2 '■ 


i-l ' 



1 4- 

^ 3 r; 



Due to the regularity of equation E2 we have 

/'lo = 0. 



(4.2.14) 



(4.2.15) 



Using Trapezoidal and Simpson's rules we arrive at the following expressions for / on 
the initial slice: 



/i = 



/o + lAz;(o-|) 

1 + iAz;fi ' 
^ '1 



/2 



/0 + lA.(0-4^-^-|) 



fr 



1 + 




















i-2 ' 


i-l 4r 



,) 



1 + lAz;| 



(4.2.16) 



(4.2.17) 



(4.2.18) 
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4.2.3 Evolution Scheme 



The evolution scheme evolves the massless scalar field coupled to Einstein's equations 
for one time-step. In section 3.8 we have chosen to use a constrained first order system 
consisting of the equations E2, SI and D2, D5, C2, D4. 

Basically, since the equations to be solved are coupled ODEs, we need to alter- 
nate prediction/ correction steps in the u- and z;-directions. We have implemented the 
integration of the u-equations via the explicit Trapezoidal Runge-Kutta method, equa- 
tion (C.2.8), while we use the implicit Trapezoidal method, equation (C.2.7), for the 
v-equations to render the scheme stable. The set of discretized equations forms a non- 
linear coupled algebraic system where some of the predicted function values need to 
be used in the correctors of other equations. The predicted value for a gridfunction / is 
denoted by the hatted quantitiy /. 

First, we tackle the boundary conditions at the center of spherical symmetry r — 0. 
They can be categorized as follows: 

• algebraic boundary conditions for r, g, and f 

• homogeneous differential boundary conditions for a and s involving first deriva- 
tives |b 

• nonhomogeneous differential boundary conditions for p=q involving first deriva- 
tives lit 

Depending on the number of past slices that have already been calculated, we could 
discretize the differential boundary conditions in a number of ways. For first derivatives 
at the origin we could use the following stencils: 



'^^ 
dr 

. ds 



-3fl+ -I- 4fli - fl, , 

= + 0{^r^) (4.2.19) 

f,r=o 2Ar ^ ' ^ ' 

3s+ - 4s,' -I- ST 

= ' 2Am ' +^(^" ) (4.2.20) 



It would, however, be wrong to use those approximations, since the grid will not, in 
general, be equidistant in r. Therefore we must resort to more general finite difference 
formulae as has already been discussed in section 4.2.1. 

We assume that indefinite expressions (at r = 0) of the form 0/0 which occur in 
the right-hand sides of some of the ODE's have already been regularized. (Again, see 
section 4.2.1.) 

The core evolution scheme uses RK2 (explicit Trapezoidal) for the evolution equa- 
tions in the u-directions (the d- and ij-equations), and implicit Trapezoidal integration 
for the c-equations. Luckily, almost all of the resulting coupled non-linear algebraic 
equations in v can trivially be made explicit by solving them in a fixed order. Only the 
g- and r-equations need to be decoupled by solving a linear system (See Ref. [Har05]). 

Effectively, we first calculate predictor values d and q (via an explicit Euler step) for 
the M-equations. Then, we proceed to calculate s, a, g and f, and finally / and p, which 
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are already the final values for those quantities. Last, we correct d and q according to 
the RIC2 explicit trapezoidal step and set the final values equal to the hatted values for 
s, a, g, r, f and p. We give a detailed description of the algorithm in pseudocode: 



Evolution Algorithm: Assuming that the boundary conditions have already been 
enforced at the origin gridpoint with index iO'^ , we loop over all remaining 
gridpoints in a w = const > slice: 

1. Predictors 



• w-equations: Euler forward 



d = di + Au\ r piqi 



q = cji- Am 



fi^i + giPi 



(4.2.21) 
(4.2.22) 



• ?7-equations: implicit Trapezoidal 



s = 



s = 



1 

s + -Avq 



a = 



1 - \Avd 



g = gU + lAv(2dl,gl,-rUql,f) 
1 



r = 



g = 



r = 



g - ^Avq^f 
l-Avd + (^f 
^Avg + (1 - Avd)f 



(4.2.23) 
(4.2.24) 
(4.2.25) 
(4.2.26) 

(4.2.27) 
(4.2.28) 

(4.2.29) 
(4.2.30) 



If we are at the first gridpoint out from the origin 



/ - fiLi + 2 ^^/regularized RHS 



(4.2.31) 



64 I Chapter 4: Numerical Algorithms 



else 



1 /^i^^i 



i-l 



/ = 



/- \^vcP■|f 
1 + \^vg|f 



(4.2.32) 
(4.2.33) 



I£ we are at the first gridpoint out from the origin 



(4.2.34) 



else 



'i-l 



V = 



p - \Avfq/f 
1 + ^Avg/f 



(4.2.35) 

(4.2.36) 
(4.2.37) 



2. Correctors 



• M-equations: explicit Trapezoidal 



d + + Aw ^ pq 



q + qi- Am 



(4.2.38) 

(4.2.39) 
(4.2.40) 



• y-equations: trivial 



sl=s 
aj - a 

81 =S 
rt - f 



f--f 

+ ^ 

Pi =P 



(4.2.41) 
(4.2.42) 
(4.2.43) 
(4.2.44) 
(4.2.45) 
(4.2.46) 
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4.2.4 Horizon Detection and Excision 



To detect apparent horizons we employ a utility function that reliably finds a possible 
zero in the gridfunction g [i] that represents g{v) = ri,{v) for a given null slice u — const. 
This is accomplished by looking for sign changes in the discrete gridfunction via a linear 
search. Such an algorithm is straightforward to implement and safe to use, as it cannot 
overlook zeroes. As we have discussed in section 2.4.1, due to our use of an outgoing 
null foliation and the causal properties of the outer trapping horizon, the function g{v) 
can, on a given null slice, at most have one zero, so that we can safely ignore any further 
zeroes as numerical artifacts. 

In the code, we compute a quantity r_AH defined as the average of r[i_left] and 
r [i_right] , the interval bracket of the root of g [i] . I.e. the actual zero vah is contained 
in the interval [r[i_left] ,r[i_right]]. 

We illustrate the behavior of the function g by some results from numerical evolu- 
tions. Figure 4.7 shows g{v) for some null slices 
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Figure 4.7: This figure shows gip) - r' for a series of m = const null-slices. The first three slices 
are untrapped, while the last three contain an apparent horizon, the radius of which decreases 
with increasing retarded time u. 

Owing to the properties of the outer trapping horizon (see section 2.2), one might 
assume that in a numerical collapse evolution using an outgoing null slicing the first 
AH is always detected at the outermost gridpoint of the slice. In practice, this is simply 
not the case. Due to the finiteness of the timestep, Au, it is highly unlikely to exactly 
"hit" the null slice Meh which coincides with the event horizon. Instead, one would 
expect to "jump" over the event horizon and detect an AH for a u > meh arid thus, since 
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0.0001 0.001 0.01 



r 

Figure 4.8: Instead of showing the null expansions d± which diverge at the origin due to their 
1/r dependence, we depict f = r and g = r' which decide the sign of the expansions as functions 
of r. We show one slice shortly before apparent horizon formation (red and blue curves) and 
a second slice containing a marginally trapped surface (green, violet). Note that the function / 
barely changes between the two slices, while g abruptly bends down into negative values. In 
contrast to figure 4.7, where we plotted g as a function of advanced time v, g{r) first grows with 
r, and after reaching the horizon radius, it decreases again; the null slice "bends" back towards 
the curvature singularity in the black hole. The condition for an apparent horizon, 6+ - while 
0_ < is fulfilled. 
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Figure 4.9: This diagram shows the location of the trapping horizon in a numerically generated 
spacetime by the double-null code. The coordinates are linear combinations of the null coor- 
dinates {u,v). This choice allows us to map null slices into diagonal lines. The trappped and 
untrapped regions have been shaded light and dark blue. The trapped region is bounded by 
the trapping horizon and the outer grid boundary. In addition, we indicate an approximation 
of the event horizon. 
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z^AH decreases with increasing u (if the AH is spacelike, which it has to be if matter falls 
in), we detect the AH somewhere inside the grid, but not at the outer grid boundary. On 
a slice that penetrates the trapping horizon, all gridpoints for which v > vah are trapped 
and the slice bends back towards the curvature singularity at r = 0, i.e. r decreases until 
it reaches r = for some v > Vah- Effectively, this forces us to cut off a sizeable chunk of 
the grid on each slice that contains trapped surfaces: We look for trapped gridpoints, 
g[i] < 0, that have fallen into the curvature singularity, that is they fulfill r[i] < 0, and 
excise them and all other gridpoints that we deem too close to the singularity. This is 
simply done, by setting the outer grid boundary index N to a new value (smaller than 
the old one) such that all offending gridpoints are removed from the grid. In contrast 
to Cauchy evolution codes, excision is very natural and straightforward to implement 
in a double-null code. 

4.2.5 Mesh Refinement 

Although we closely follow the work of Hamade and Stewart [HS96] in our imple- 
mentation of the double-null code, we employ Garfinkle's [ ''^^] approach to increase 
resolution which we already described for the DICE code in section 4.1.3, rather than 
the more involved Berger Oliger type adaptive mesh refinement [B084], which was also 
used in the seminal work by Choptuik [d'192a, Cho93]. The application of Berger Oliger 
type AMR in a characteristic framework has been discussed by Pretorius and Lehner 
[PL04]. 

Concerning the mesh refinement strategy, the only difference to the DICE code, is 
that the double-null code is uncompactified and we can adjust the outer boundary of 
the grid to be slightly larger than the null ray v - v* which hits the accumulation point 
at {u*,v*) of the self-similar solution as described in figure 4.10. Self-similar solutions 
and critical collapse phenomena are discussed in detail in chapter 5. 

4.2.6 Constraints and Convergence 

To check that our numerical solution fulfills the Einstein equations we calculate the 
constraint CI: 

/ = 2^/ - rp2 (4.2.47) 

and suitably normalize it with the absolute value of the running maximum of its indi- 
vidual terms over a slice, as described in section 4.1.4. We also compute the discrete 
f2-norm of the normalized constraint. Figures 4.11 and 4.12 show the absolute value 
of the constraint CI as a function of r for a couple of null slices in a subcritical and a 
supercritical evolution, respectively. 

Along with the monitoring of constraint equations, convergence testing is an im- 
portant strategy to ensure the correctness of numerical results, namely that, as the 
gridspacing tends to zero, the numerical approximation converges to the analytical 
evolution system according to the accuracy of the discretization. The basic methodol- 
ogy is explained in appendix D. 
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Figure 4.10: This figure shows how to adjust the outer grid boundary fouter (or equivalently 
^'outer) in order to accurately resolve the self-similar features of critical collapse close to the origin. 

In the following, we present results from a convergence test using 1001 gridpoints as 
the base resolution and 2001 and 4001 gridpoints for the medium and finest resolution, 
respectively. As the initial data are not as close to criticality as those used in the 
convergence test shown in figure 4.2 for the DICE code, it was not necessary to adjust 
the initial data parameter p for each resolution. Due to the size of the gridspacing and 
second order accuracy of the code, the finite difference error is approximately 10~^ for 
the coarsest, 2 x 10"'^ for the medium and 6 x 10~^ for the finest resolution. The set of 
evolutions under consideration forms an apparent horizon at approximately ii = 1.088 
with a mass of 0.0076 + 0.0001. For details see table 4.1. 



N 


uah mAH 


1001 
2001 
4001 


1.088 0.00775767 
1.088 0.00765987 
1.088 0.00763467 



Table 4.1: The initial AH masses (i.e. when we first detect an AH), for the three evolutions of 
the convergence test. 



As can be seen from figure 4.13 the AH is well resolved even for the coarsest 
resolution of 1001 gridpoints. A convergence test for the /'2-norm of the constraint CI 
is shown in figure 4.14. Interestingly, the first part of the evolution exhibits fourth 
order convergence. Later on we have solid second order convergence. The constraint 
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Figure 4.11: Here we show the spatial variation of the absolute value of the constraint CI for 
some slices of a subcritical evolution. The maximum of |C1| over the course of the evolution is 
below 10~^. Figure 4.12 shows the constraint for the supercritical case. 
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Figure 4.12: This figure shows the absolute value of the constraint CI for a few slices of a 
supercritical evolution. On the null slice # 2175, an apparent horizon has formed at r ^ 0.015 
which causes the constraint to grow about three orders of magnitude. 
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Figure 4.13: This figure shows the mass-function m{r) for two successive null slices in the lowest 
resolution (1001 gridpoints) run of the convergence test. An apparent horizon is detected on 
slice # 544 at u = 1.088. The detail shows that the AH is well resolved. 
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Figure 4.14: This figure shows the convergence of the /'2-norm of the constraint CI. Initially, 
the convergence is fourth order, but note that the finite difference error for the lowest resolution 
run IS about 10"^ For I/Mbh > 100 convergence is roughly second order. An apparent horizon 
formes at m = 1.088 with a mass of 0.0076 ± 0.0001. The constraint violations reach 0.1 at the 
formation of the apparent horizon, but otherwise are fine. 

violations reach 0.1 at the formation of the AH, which is the most demanding part of a 
collapse evolution, but otherwise are fine. In figure 4.15 we show that the scalar field s 
is second order convergent. 

4.2.7 Timelike Observers 

In the following, we describe how to calculate the proper time of an r = const observer 
from the expression (3.8.35) derived in section 3.8.3. As the r-values of our gridpoints 
on the double-null grid keep changing and we need to integrate over time, we fix a set 
of predetermined radii r [i] and use interpolation to obtain the corresponding set of v 
coordinates for each u = const slice. We store the running integrals in a vector tp[i]. 
The vector is initialized at the initial slice by setting tp[i] = r[i], since we assume a 
flat geometry, initially. On each null slice, and for each observer at r = r [i] we 

1. Determine and 5 via the relations 



r{u,v) = r[i] 



(4.2.48) 




(4.2.49) 
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2. Evaluate 



dv 

du 



=r[i] 



V - V 

Au 



(4.2.50) 



to first order accuracy. 

3. Evaluate a{u, v) using spline interpolation. 

4. Update the running integral by adding the term Au a{u, v) 




u — Am 



r = const 



Figure 4.16: Calculation of |^|^ 

If the slice contains an apparent horizon, things are more complicated, as t{u, r) is no 
longer a single-valued function of r. We first compute the index of the gridpoint with 
the maximum r-value in the grid. We can then use spline interpolation in each of the 
two halves of the grid, where r{u = const, v) is either an increasing (untrapped region) or 
decreasing (trapped region) monotonic function of v or gridpoint index. In the trapped 
part of a slice ^|^. < 0, so we take its absolute value in (3.8.35). We now have to keep 
track of one running integral for proper time for each half of the slice. 

For near-critical evolutions, the trapped part of the slices is very small and hardly 
warrants the difficulties incurred. Moreover, due to excision of gridpoints that come 
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close to the singularity at r = 0, we are missing data to compute tp [i] in the outer half 
of the slice. 
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Chapter ^ 

Critical Phenomena 



In general relativity, critical phenomena can be subsumed by three essential features: 
scaling, self-similarity and universality. They arise in the behavior of near-critical 
evolutions, in the course of which a universal self-similar solution is approached for 
some time and the idiosyncrasies of the initial data are mostly "forgotten" and replaced 
by symmetries inherited from the attractor. Ultimately, the self-similar solution causes 
the evolution to approach an endstate of the system. While self-similarity can be 
analyzed in single evolutions, scaling only becomes apparent when considering the 
final parameters of a one-parameter family of evolutions. Dimensionful quantities, 
such as the black hole mass, exhibit a power-law behavior with an exponent that is 
related to the unstable mode of the self-similar solution. 

We first define the special symmetries, discrete and continuous self-similarity, that 
lie at the core of critical phenomena. Next, we give a qualitative description of near- 
critical evolutions in the dynamical systems picture. Finally, we derive the scaling law 
for dimensionful quantities. 



5.1 Self-Similarity in General Relativity 

First, we introduce the concepts of discrete (DSS) and continuous (CSS) self-similarity 
and then go on to talk about coordinates adapted to these symmetries. Although we are 
not interested in CSS per se, it is often much easier to handle and useful for modelling 
situations where the periodicity of DSS does not play a decisive role. 

A spacetime (M,g) is said to be discretely self-similar (DSS) [Gun99] if it admits a 
discrete diffeomorphism : M — > M which leaves the metric invariant up to a constant 
scale factor: 
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where A is a dimensionless real constant, n £ N and 



(5.1.2) 



is the pullback of the spacetime metric under the diffeomorphism Oa and we denote the 

Jacobian of this map by 

Similarly, a spacetime that admits a one-parameter family of such diffeomorphisms, 
parametrized by A, with Oq being the identity, is called continuously self-similar (CSS). 
The generating vector field = ^^a|^_q is homothetic, that is, it obeys the conformal 
Killing equation with a constant coefficient on the right hand side: 

(5-1-3) 

The choice of the constant is pure convention. 

Following Gundlach [Gun97b] we introduce coordinates (t, x'^) such that if a point p 
has coordinates {t,x"), its image Oa(p) has coordinates (t - A,x"). One can then verify 
that DSS in these coordinates is equivalent to the condition 

gab{T,x") = e^-'g,,{T,x''), (5.1.4) 

where 

(5.1.5) 

8abiT,x")^gab{T + A,x''). (5.1.6) 

To explain the connection between CSS and DSS one may introduce a vector field 
^ = d/dj. The discrete diffeomorphism is then realized as Lie-dragging along E, by 
a distance A. The discrete diffeomorphism cannot uniquely determine ^, as, loosely 
speaking, the vector field f is free to do whatever it wants between DSS "echoing" 
periods. 

We introduce adapted coordinates (t, z, 6, (p) based on the Bondi coordinates defined 
in section 2.3.2: 

T = -In^^-^, (5.1.7) 

T re^ 

(5.1.8) 



where u* is a real number which denotes the accumulation time of DSS and C(t -I- A) = C,{t). 
A detailed derivation of these adapted coordinates is given in [LecOl]. We also require 
that C,{t) > C(t) for all t such that ?'t(t,z) < for all t. The remaining gauge freedom in 
C(t) is chosen such that dj (which is timelike at the origin) becomes null at z = 1. This 
null hypersurface, the backwards lightcone of the accumulation point, is known as the 
past self-similarity horizon (SSH). See figure 5.1 for a spacetime diagram that illustrates 
the behavior of these coordinates. 
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Figure 5.1: This diagram shows a prototjrpical self-similar spacetime. The adapted coordinates 
(t,z) defined in equations (5.1.7) and (5.1.8) cover only the region u < u*. The metric functions 
/3, J and m are constant along lines of constant z. Thus, their features shrink to zero as one 
approaches the accumulation point at (m = u*,r - 0). Self-similarity is confined to the backwards 
lightcone of the accumulation point, which extends from the origin to the past self similarity 
horizon at z = 1 . 
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Next, we transform the Bondi line-element (2.3.17) to these new coordinates using 
the inverse transformation 



which yields 



^(t,z) + 2z(C(t)-C(t)) 

- 2e^f^^^'^^u^Q{T)dTdz + [uzQ{z)f dQ^j. (5.1.11) 

This line-element fulfills the condition (5.1.4) if the metric functions |S and V/r are 
periodic in t with period A: 



u = u*(l- e-^) 
r - u*ze~'^C,{T), 



dx" 



(5.1.9) 
(5.1.10) 



ig(T,z) = |g(T + A,z) 

V V 

-{t,z) = y(T + A,z). 



(5.1.12) 
(5.1.13) 



Which conditions does the symmetry (5.1.1) imply for the scalar field? Clearly 
matter has to share at least some of the symmetry of the geometry as it is coupled via 
Einstein's equations. First, we need to know how the Einstein tensor behaves under the 
action of the discrete diffeomorphism. Thanks to general covariance, the pullback of 
the Einstein tensor is equal to the Einstein tensor formed from the pulled-back metric: 



o;G[^] = G[o;g]. 



(5.1.14) 
Because 



It follows from the definition of DSS, equation (5.1.1), that O^G[^] = G ^^'^g 
the structure of the Einstein tensor is such that in every term that contains the metric, 
the inverse metric also appears, and since ^\g~\ - ^'"^^glp we find G \^*/^g\ = G[g] and 
thus 

%G[g] = G[g]. (5.1.15) 
Via Einstein's equations this symmetry also holds for the stress energy tensor: 

Taking the pullback of the scalar field stress energy tensor, equation (3.1.4), yields 



«fclp dx" dx^ 
dx" dx'' 



1 d(t>'^ d(t>'^ 



(5.1.17) 
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where we have used the pushforward of the inverse metric. 



and the definition of DSS in terms of the metric, equation (5.1.1). To simplify the 
expression we use the adapted coordinates (5.1.7) and (5.1.8), so that the point p has 
coordinates (t, z) and Oa(p) has coordinates (t - A, z) and drop the angular dependence: 

V,(/)(T,z)Vi,(/)(T,z) - ^^„fo(T,z)/^(T,z)V,(/)(T,z)Vrf(/)(T,z) 

= V«(/)(t - A,z)Vi,(/)(T - A,z) - ^gabiT,z)g'^{T,z)Vc(p{T - A,z)V^(/)(t - A,z). (5.1.19) 

Therefore, we find the condition 

V„(/)(t,z) = +V,(/)(t-A,z). (5.1.20) 

For the scalar field, the most general ansatz compatible with DSS is (see Ref. [Gun97b]) 

(f){x') = c(){t,z) + KT with (p{T + A,z) = +(^(t,z). (5.1.21) 

The constant k is determined by the (uu) component of Einstein's equations written in 
terms of the adapted coordinates (t,z). As discovered numerically by Choptuik, k = 0, 
for unknown reasons. 

The vanishing of k together with the choice of the minus sign in equation (5.1.21) 
for cf) gives rise to a further symmetry (see [Gun97b], [LecOl]) that has been numerically 
observed: If the metric functions are periodic with a period A and the scalar field is 
"anti-periodic" with respect to this period, then, the field is obviously periodic with 
respect to the period 2A. We say that the solution is DSS with period A = 2A and has 
the additional symmetry 

/*(T + nA/2,z)=f(T,z), (5.1.22) 

where /* denotes the metric functions |S, V/r, m/r and the function C,{t), while the scalar 
field (p satisfies 

KJ'^ = (-!)>• (5-1-23) 
This is consistent with the fact that |S, V/r, m/r are even in (p. 



5.2 General Relativity as an Infinite Dimensional Dynamical 
System 

Following [Gun99], we give a brief qualitative introduction to the dynamical systems 
picture. We pretend that general relativity can be treated as an infinite-dimensional 
dynamical system and ignore problematic issues such as convergence to an attractor. 
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We restrict ourselves to isolated self-gravitating systems, such as a ball of radiation or 
a star. The phase space consists of the space of initial data for an isolated ^ system, which 
in our case is just the space of (p{uQ,r). Possible end-states of the system correspond 
to basins of attraction in phase space. The boundary between two basins of attraction is 
called a critical surface. We restrict attention to the existence of two distinct endstates; 
otherwise, more complex critical points would arise at the intersection of more than two 
basins of attraction. Minkowski space is an attractive fixed point in the dispersion basin, 
while Schwarzschild black holes form a half-line of attracting fixed points, parametrized 
by mass. 

Consider an attracting fixed point or limit cycle, the critical solution, that lies in this 
critical surface, as shown in figure 5.2. We assume that this critical solution, in turn, is 
an attractor of codimension one in a neighborhood of phase space, which means that 
it possesses a single growing perturbation mode transverse to the critical surface and 
an infinite number of decaying modes tangential to the surface. ^ Trajectories starting 
out in the vicinity of the critical surface at first move close to the surface towards the 
critical solution. In the so-called intermediate asymptotics (a term coined by P. Bizon) they 
spend some time in the vicinity of the universal critical solution, shedding all the details 
of initial data from which they originated. After a while the growing mode becomes 
dominant and the trajectories depart towards the flat space or black hole fixed point 
depending on whether the data were sub- or supercritical, respectively. 

In the context of this work the terms "critical solution" and "DSS solution" (or 
"self-similar solution") are often used interchangeably. Due to self-similarity the DSS 
solution is evidently not asymptotically flat. It is however not so clear whether the term 
critical solution entails asymptotic flatness or not. Here, we want to make this sloppy 
usage more precise: We mean by critical solution an asymptotically flat spacetime which 
agrees in a finite region near the centre with the DSS solution, but falls off at null infinity. 
E.g. it would be consistent to have the exact DSS spacetime in the region defined within 
the backward light cone of the DSS singularity (the past self similarity horizon) and 
outside to fall off smoothly to infinity. Similarly, we call a solution near-critical, if it 
comes close to the critical solution (in the above sense) near the centre of collapse for 
some time during the evolution. The spacetime region in which an evolution spends 
in the neighbourhood of the critical solution depends on the fine tuning and the class 
of initial data. Therefore, the failure of the DSS solution to be asymptotically flat is 
perfectly consistent with its role in the dynamics of a localized object emitting radiation 
to null infinity. 



^ Note that DSS or CSS do not lie in this phase space since they are not asymptotically flat. Solutions 
which are only self-similar inside a finite radius and asymptotically flat can be included in the phase space. 

^ Attractors which possess more than one imstable mode would be exceedingly hard to "run into" in 
the course of numerical evolutions, as they would require the concurrent fine-tuning to "critical values" 
of more than one parameter. 
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Figure 5.2: A qualitative picture of the phase-space of near critical evolutions. Every point 
corresponds to one configuration (p{u = const, r). The stable manifold contains a codimension 
one discretely self -similar solution, depicted as a limit cycle, as it is periodic in t. We have 
indicated a one-parameter family of initial data by a line that intersects the critical surface on 
the left. At the "critical" parameter value p = p* , the trajectory is attracted by the stable modes of 
the critical solution and never leaves the stable manifold, ultimately arriving at the DSS solution 
in the limit t ^ oo. Trajectories with p > p* may approach the critical solution for some time, 
as long as the admixture of the unstable mode is small, but ultimately leave the vicinity of 
the critical solution and move on to form a black hole. In contrast, configurations with p < p* 
eventually disperse and reach flat space. 
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5.3 The Black Hole Mass-Scaling Law 



In this section we present an argument for the scaling of the black hole mass. We follow 
Gundlach ([Gun99], [Gun97h]) and Lechner [LecOl]. 

For the sake of simplicity, we first consider the spherically symmetric CSS case. 
Assume that the self-similar solution is an attractor of codimension one in phase space, 
that is, it has exactly one unstable mode with eigenvalue Aq. Furthermore, assume that 
the system has two possible endstates, dispersion to flat space or black hole forma- 
tion, such that the stable manifold of the self-similar solution divides phase space into 
subcritical and supercritical data. (See figure 5.2 for an illustration.) Let Z denote any 
scaling variable (j3, j,(p, 7). We use adapted coordinates to CSS: 

T = -]n- and z- — - — . (5.3.1) 

M* M* - U 

General near-critical initial data are attracted via the stable modes of the CSS solution, 
which we denote by Z,(z). In contrast to a DSS solution, which would be periodic in the 
adapted coordinate t, a CSS solution is independent of t. In the echoing region, where 
Z» dominates we write the solution as a small perturbation of the CSS solution: 

Z(t,2) = Z,(2) + Co(p) e^^' Zo(z) + 6Zstable(T,z). (5.3.2) 

We assume that the single unstable eigenvalue Aq is real and positive. As t — > oo 
all other perturbations, which are contained in 6Zstabie(T,z), vanish. The amplitude 
of the unstable mode Zo(z) contains information on the initial data, in particular, the 
parameter p that characterizes a one-parameter family of initial data. For so-called 
critical initial data, p — p* , the unstable mode is completely tuned out, so that Co(p*) = 0. 
For near-critical data, we may then linearize around p*: 



{p-f)+0{{p-ff), (5.3.3) 



so that for suitable t 



Z(t,2)-Z,(z) + ^ 



{p-f)e'°'Zo{z). (5.3.4) 



Clearly, as p — > p*, we see more and more of the critical solution. Equation (5.3.4) also 
sheds light on the universality of the self-similar solution: The critical phenomena are 
independent of the family of initial data. 

Let Tp be a value of the adapted time chosen such that the stable modes are negligible 
compared to the unstable mode and that the amplitude of the uristable mode equals a 
small number e «c 1: 

^ {p-p')e^°^p =e. (5.3.5) 
dp p. 

Note that Zp depends on the initial data parameter p. We rewrite this relation as 

e-'p = const {p - p*f'^° , (5.3.6) 
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where the constant is independent of p. 

For large t the linear approximation eventually breaks down when the unstable 
mode has grown beyond a certain bound. If, as we assume in the following, the data are 
supercritical, a black hole forms at late times. We emphasize that we need not follow 
the evolution into the nonlinear regime; it is sufficient to note that intermediate data 
which we can extract at t = Zp depend on r only through the coordinate z: 

Z(Tp, z) ^ Z,(z) + eZo(z). (5.3.7) 

In Bondi coordinates we have 



Z(wp, r)^ZJ + eZo . (5.3.8) 

The crucial point of this argument is that these intermediate data at u - Up depend 
on the initial data, say at time u — 0, only through the overall scale 

L^u -Up = ue^'p. (5.3.9) 

The field equations do not have an intrinsic scale so that in the absence of an external 
scale, such as a cosmological constant, the field equations are invariant under rescalings 
u — > rju, r — > 7]r. Therefore, the solution based on the data at u = Up must be universal 
up to the overall scale. 

Because the black hole mass has dimension of length in geometrized units (G = c = 
1), it must be proportional to L: 



Mbh oc L cc (p - p*)l/'^° 

and we define the critical exponent 



7 := 



Re Ao 



(5.3.10) 



(5.3.11) 



Relation (5.3.10) was derived independently by Koike et al. [KHA95] and Maison 
[Mai9b]. 

As pointed out by Garfinkle and Duncan [ 7D98], it is also possible to analyze 
scaling in subcritical evolutions via the scaling behavior of the maximum over a whole 
evolution of the Ricci scalar curvature at the axis. Due to its dimension being 1/L^, the 
Ricci scalar scales as 

maxR(M,0) ~ (p - p*)-^'^° . (5.3.12) 

Numerically, it is easier to measure the critical exponent from the subcritical scaling 
of the Ricci scalar than from the supercritical scaling of the black hole mass, since one 
avoids numerical difficulties incurred in apparent horizon formation. 

The mass-scaling relation (5.3.10) is universal in the sense that it holds for all one- 
parameter families of initial data, since its derivation relies on perturbations of the 
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universal (for each model) critical solution. Thus, all solutions which approach the 
attractor for a while will exhibit this same scaling law. 

Although the critical solution itself is not asymptotically flat due to its self-similarity 
one can restrict the analysis to a region of a sufficiently large, but finite radius, so that 
the black hole mass is not affected. 



5.3.1 DSS Fine Structure 

In the case of a discretely self-similar critical solution, the scaling law (5.3.10) is modified 
by an oscilliatory "fine structure" overlaid on the linear relation between InMsH and 
]n{p — p*). This phenomenon was first predicted by Gundlach [Gun97b] and Hod and 
Piran [HP97]. 

Similar to the CSS case, where the perturbation modes of the background solution 
would retain its CSS symmetry (independence of t), the modes are now periodic in t 
and equation (5.3.2) becomes 

Z(t,z) = Z,(t,z) + Co(p)e'^°"Zo(T,z) +6Zstable(T,z). (5.3.13) 

Note that 6Zstabie(T^/2) = LSr=i Cne'^"^Z„(T,z) is not periodic in t, whereas the individual 
modes Z„ are, since the perturbation must break the DSS symmetry or the critical solu- 
tion would not be isolated [Gun97b]. Again, for suitably large t so that the perturbation 
is still small and we can neglect the stable modes, we obtain for near-critical data 



Z(t,z) -Z,(t,z) + 



dCp 
dp 



{p-p*) e'^o^Zo(T,z) 



(5.3.14) 



As before, we extract intermediate data at a time defined as 



Ao 



In 



(P - Pl 



e dp 



(5.3.15) 



where e is a small constant. These data depend on r only through the dimensionless 
coordinate z: 

Z(Tp,z) ^ Z,(Tp,z) -I- eZo(Tp,z), (5.3.16) 
so that in Bondi coordinates we have 



Z{up,r) ^ Zo 



(w» - Wp) C(Wp)^ 



+ eZo 



Up, 



{u* - llp)Q{tlp)^ 



(5.3.17) 



where Up — u* (1 - e ^p). The entire solution evolved from these data scales with 

(u* - Up) C(Mp) f{up) = u^e-^'QiTp) /(Tp) = e-'p /(zp), (5.3.18) 

where /, / (and Q are periodic functions of their argument. Therefore, the black hole 
mass obeys the following scaling law 



niBH = ci{p- pf^° / [- In ci - 1/Ao In (p - p)] 



(5.3.19) 
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where the constant Ci 



idCol 
e dp \f 



1/Ao 



depends on the family of initial data. Equivalently, 



In Wbh = In Ci + y In (p - p*) + /[- In Ci - 7 In (p - p*)] , 



(5.3.20) 



where / = In / and we have introduced the scaling exponent 7 = I/Aq. As a function of 
ln(p - p*), / is periodic with period A/(2y) « 4.61 (for the massless scalar field). This is 
due to the additional symmetry of the metric functions mentioned in equation (5.1.22) 
(see Ref. [Gun97b]). 
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Quasinormal Modes and Tails 



6.1 Introduction 

In general relativity, quasinormal modes (QNM) (see [KS99, FN98, Cha83] for a review) 
arise as perturbations of stellar or black hole spacetimes. Due to energy loss caused 
by radiation, these infinite systems do not exhibit normal mode oscillations which 
are characteristic of compact linear oscillating systems without damping. Instead, 
the frequencies become "quasi-normal", i.e. complex, with the real part representing 
the frequency of the oscillation and the imaginary part representing the damping. 
Moreover, in contrast to normal modes, quasinormal modes do not form a complete set 
of basis functions. 

In addition, for late times, the weak decay of the Schwarzschild potential causes the 
quasinormal oscillations to be swamped by the radiative tail[ ], which only decays 
polynomially in time. 

In the following we give a brief review of the evolution of weak fields on a 
Schwarzschild background (or similar) [NF89]. The general procedure is to expand 
the field in spherical harmonics (according to its spin s). For each radiative multipole 
/ > s there is a scalar field which depends on the Schwarzschild coordinates t and r. 

(s) 

Each such scalar function (p^ (t, r) satisfies an equation 



= Vy-'irW^' (6.1.1) 



where is the effective potential and r, is the usual Schwarzschild "tortoise" coordi- 
nate, which is related to r by dr/dn - 1 - 2M/r. 

(s) (s) ' f 

With the ansatz cp^ (?', = i/'J {r*)e~ we arrive at the following second order ODE 
that is similar to the one-dimensional Schrodinger equation for a particle encountering 
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a potential barrier on the inifinite line[Iye87]: 



+ 



(6.1.2) 



For scalar fields the potential is 




(6.1.3) 



The effective potential behaves as a potential barrier, as it is appreciably distinct from 
zero only in the neighborhood of r» « (r « 3M). As r, — > +oo the potential falls off very 
rapidly. 

The so-called quasinormal modes of the black hole are solutions to the perturbation 
equation (6.1.2) which satisfy radiation boundary conditions for purely outgoing waves 
at (spatial) infinity and purely ingoing waves at the horizon. The real part of the QNM 
represents the oscillation frequency, while the imaginary part represents the decay. 

The radiation produced in response to a perturbation of the field around a black 
hole can be divided into three stages: 

• radiation emitted directly by the source of perturbation 

• radiation due to the damped oscillations of quasinormal modes excited by the 
perturbation source ("ringing radiation") 

• power law tails of radiation, caused by scattering of waves by the effective poten- 



A distant observer first records the radiation from the perturbation source, then the 
quasinormal ringing, which decay exponentially, and finally the radiation tails which 
decay much slower, by a power law, and have a very small amplitude in comparison to 
the first two components. 

As an example, the dimensionless frequency Mcu for the I = 2,n — mode of 
gravitational perturbations is Mcu w 0.37 - 0.089f, while the / = 0, n = mode for scalar 
perturbations is Mcu « 0.11 - O.llz [Iye87]. Both are fundamental modes. 

To convert from the dimensionless frequency Mco in geometrized units (G = c = 1) 
to non-geometrized units, we note that the dimension of Mco is mass/time which yields 
the conversion factor c^/G. The QNM frequency for a 10 solar mass black hole is then 



tial 



/ 



1 Re Mco 
2nG WMq 



Hz « 3.2Re (Mcu)kHz 



(6.1.4) 



and the damping time is 




G IOMq 
c3 |Im Mco\ 



s 



4.97 X 10-2 
|Im Mcu| 



ms. 



(6.1.5) 
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For the / = 0, n = scalar mode, this corresponds to a frequency of 357Hz and a 
damping time of 0.45ms, while for the I - 2,n — gravitational mode, this corresponds 
to a frequency of 1 .2kHz and a damping time of 0.55ms. Compared to the Schwarzschild 
radius for the 10 solar mass black hole, which is rs = ^^^^ ~ 30km, the wavelengths of 
the fundamental QNM are much bigger, 841km and 250km for the / = 0, n = and the 
/ = 2, n = mode, respectively. 



6.2 Quasinormal Modes In Critical Collapse 

As we have mentioned, QNM excitations are, in general, obtained from linear pertur- 
bations off a fixed background, together with their associated (complex) eigenvalues. 
Thus, in a highly dynamical setting, such as in critical collapse evolutions, one would 
not expect to see (identify) quasinormal modes. However, it turns out that for our 
setting, the least damped spherically symmetric mode for scalar perturbations of a 
Schwarzschild black hole plays a relevant role. 

Perturbation theory[iyeb/ ] gives the following value for the half-period 

where Re cy = 0.11 is the real part of the n = 0, / = QNM. This mode has previously 
been detected in supercritical evolutions (far away from criticality) for a self-gravitating 
massless scalar field by Gundlach et al. [ " PP94b]. 

In the following, we analyze radiation signals for near-critical evolutions, where the 
notion of a fixed background mass does no longer apply. We find that the monopole 
moment of the scalar field c(mb) shows a damped oscillation with exponentially increas- 
ing frequency (see figure 6.1). Moreover, the sizes of the half-periods measured from 
one extremum to the next in c{u^) roughly agree with the half-periods obtained from the 
least damped quasi-normal mode (QNM) of a Schwarzschild black hole with a strongly 
changing "background" mass Mi,g{u^) as shown in figure 6.2. M},giih) is obtained by 
evaluating m^{u^) at the mean value between the extrema of c{u^) (which are inflection 
points of mB(MB))- 



6.3 Quasinormal Modes for the Vaidya Metric 
6.3.1 The Vaidya Metric 

The Vaidya metric[Poi04, GS04, ACS06] is a generalization of the Schwarzschild solution 
with a changing mass. 

We express the Schwarzschild metric in the hybrid outgoing Eddington-Finkelstein 
coordinates {U, r) 

ds^ ^ -(l-—)du^ - Idudr + r^dO^, (6.3.1) 
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Figure 6.1: This figure shows the scalar field monopole moment c(Mb) for a near-critical (but 
barely supercritical) evolution. The half-periods measured from one extremum to the next 
roughly agree with the prediction of perturbation theory shown in figure 6.2. 
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0.001 



Figure 6.2: The exponential decay of the QNM half-periods predicted by perturbation theory 
is shown with annotated values at the midpoints between the points of inflection for the same 
near-critical evolution as in figure 6.1. 
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where the null coordinate u = t - is defined in terms of the tortoise coordinate r» = 
r + 2Mln(r/2M - 1). Then we allow the mass parameter M to become a function of 
retarded time: M m{u). The resulting metric is the outgoing Vaidya metric 

ds^ " ~ ~ J ^ -2 _ 2^u^r + r^dCl^. (6.3.2) 

The only nonvanishing component of the Einstein tensor is Gaa - -(2/r^) (dm/du). For 
the metric (6.3.2) to be a solution of the Einstein equations the stress-energy tensor must 
be of the form 

where /« = -Vfltl is a radial null vector. This stress-energy tensor describes null dust, 
a pressureless fluid with energy density p = -\ I {Anr^){dm I du) moving with a four- 
velocity All the standard energy conditions [Poi04] are satisfied by Tab if dm/du < 0. 
This solution of Einstein's equations describes a unidirectional radial flow of unpolar- 
ized radiation in the geometric optics (high frequency) approximation. 
In contrast, the ingoing Vaidya metric is given by 

ds^ = - (l - J dv^ + 2dvdr + r^dOi^. (6.3.4) 

and m must be a monotonically increasing function of advanced time v — t + r^^ in order 
to satisfy the energy conditions. 

Figure 6.3 shows an imploding shell of massless radiation (modelled by the ingoing 
Vaidya solution) which forms a black hole. For a given value of retarded time, u, let the 
shell be contained in the interval [t'O/^'i] of advanced time, i.e. let m{v) = for v < Vq 
and m{v) - M for v >Vi. Assume that in the interval [vq, Vi] the mass-function is chosen 
such that no shell focussing singularity appears [Kur84]. ^ For v < Vq observers are 
unaware of the infalling shell of radiation. For late times u observers close to r = will 
experience tidal stresses and will begin to fall inwards. By that time, they are already 
trapped in the event horizon. The apparent horizon of the Vaidya spacetime is always 
located at r = 2m{v) [Poi04]. 



6.3.2 QNM Modes in a Vaidya Background 

The analysis of QNM for time-dependent situations, such as accretion processes or 
Hawking radiation is discussed for the Vaidya metric by Abdalla, Chirenti, and Saa 
in [ v_sUo]. Building on earlier work by Girotto and Saa[^~' ' ], they employ a semi- 
analytical approach in which one considers the Vaidya metric in double-null coordinates 
{u,v) ah initio. After choosing a smooth mass-function m{v), it is possible to numerically 
reconstruct the metric functions. Then, scalar (or electromagnetic) perturbations on the 

^According to Kuroda[Kur84], possible choices of mass-function that do not generate a shell focussing 
singularity are m{v) ~ j^v", where w < 1, or n = 1 and j.i > 1/16. 
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Figure 6.3: This figure shows a Penrose diagram of null dust collapse. For v < vq the spacetime 
is flat. In the interval Vq < v < Vi the ingoing massless radiation described by the Vaidya metric 
forms a black hole. For v > v\ the inflow of matter has stopped and this part of spacetime is 
isometric to the exterior region of Schwarzschild. 

constructed Vaidya background are evolved using a characteristic algorithm. Analysis 
of the perturbation field close to the event horizon then yields estimates for the real and 
imaginary parts of the QNM. 

The authors point out the existence of a stationary adiabatic regime where the real part 
of the QNM varies inversely with the mass-function, just as one would obtain when 
modelling the Vaidya solution as a series of Schwarzschild slices with changing mass 
as we have done in section 6.2. In addition, they formulate a heuristic criterion that 
indicates when to expect the appearance of non-stationary behavior: 



The non-stationary behavior manifests itself as an inertial effect in the real part of the 
QNM: Re a){v) no longer behaves as \lm{v), as one would expect for a stationary 
adiabatic regime. 

In the following, we apply their criterion for non-stationary behavior to our nu- 
merical results. Note that we measure the QNM at and we are dealing with a 
monotonically decreasing function, the Bondi mass Mb(mb). ~ In section 6.2 we have 
observed that the real parts of the fundamental scalar field QNM roughly agrees with its 
Schwarzschild value (for a decreasing mass) during critical collapse. We have verified 
this (see figures 6.1, 6.2) by comparing the half-periods measured from one extremum to 

^The DICE code from which these results have been calculated cannot penetrate apparent horizons, 
and, at least for critical collapse, we have not been able to read off QNM close to the event horizon. 



\m"{v)\ > \\moj{v)\ 



(6.3.5) 
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the next in c{ub) with the half-periods predicted by perturbation theory with a strongly 
changing background mass Mijg{uB). 

We have also tried to read off the imaginary parts of the QNM from c{ub), but 
to no avail: One the one hand, this is very difficult, because we do not have many 
oscillations to accurately measure the decay. But, more importantly, as shown in figure 
6.4, c(t) behaves as This is due to self-similarity, since the monopole moment c has 
dimension of length. We note that we have corrected a slight offset in c{ub). 
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Figure 6.4: The monopole moment of the scalar field decays in accordance with self-similar 
scaling. We have corrected a small offset in c and plot |c(t) - 0.003383|. 



Instead of directly reading off the imaginary part of the QNM, we can still estimate 
the decay from perturbation theory for a few background masses and compare the result 
to d^M/du^ in order to check that the non-stationarity criterion mentioned above, (6.3.5) 
is not fulfilled. As can be seen in the figures 6.5 and 6.6, stationarity is well satisfied 
with the exception of the very late stages of collapse, where we leave the self-similar 
regime and approach black hole formation. 



98 I Chapter 6: Quasinormal Modes and Tails 



• • ln|Ima;| 




-10 t 



Figure 6.5: This figure checks that the non-stationarity criterion for QNM for critical collapse 
simulated by the DICE code. The connected filled circles represent In \coi\ at the extrema of c{ub), 
while the lower curve shows In |d^M/dMg|. With the exception of the very late stages of collapse 
we are in the stationary regime. Figure 6.6 shows a detail of this plot. 
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Figure 6.6: This figure shows a detail of 6.5 
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6.4 Power-law Tails 



As has been established by Price [Pri72], perturbation fields outside of Schwarzschild 
black holes die off with an inverse power-law tail at late times. In contrast to quasinormal 
modes this behavior does not depend on the details of the collapse process, but only 
on the asymptotic falloff of the effective potential; (i.e. in a curved spacetime, wave 
propagation is not confined to the light cones, rather waves spread inside the light 
cones, due to scattering off spacetime curvature.) Therefore, tail phenomena can be 
observed independently of the endstate of the evolution. For Gaussian initial data the 
Newman-Penrose constant [GWS94] vanishes and for late times perturbation theory 
[GPP94a, GPP94b] predicts that the field falls off as (/) oc u~-^ near J^"*" and (p oc near 
timelike infinity 

To be more precise [Bar99a, Bar99b], consider a distant static observer at r = const 
and Bondi time Amb elapsed since the "main pulse" of radiation has reached the observer 
(the duration of the main pulse has been assumed negligible in [Bar99a, Bar99b]). Null 
infinity is then found to be approximated by the region where Aug «c r within the 
context of a perturbative analysis of tail behavior [Lea86a, Lea86b, Bar99a]. This regime 
has been termed the "astrophysical zone" by Leaver [Lea86a, Lea86b, Bar99a]. 

In addition, the convergence of the perturbation expansion in [Bar99a, Bar99b] re- 
quires Amb » M, where M is the mass of the background. In our case, the mass which 
gives rise to the effective potential is bounded from above by the mass of the initial 
(ingoing) pulse, M,. Therefore, we demand that Amb » M,. Numerically, we observe 
tails only for Amb > W^Mj. 

Closeness to timelike infinity, on the other hand, demands Aub » where r* is 
the usual "tortoise" coordinate - r + 2Mln(^). In figure 6.7 we show power- 
law exponents determined by fits of xp at different r = const curves over a series of 
time intervals for a subcritical evolution. As described in section 4.1 we use spline 
interpolation to calculate ipix - const). The exponents have been determined by fitting 
the field \p at x = constant against a power of Wb in 5 distinct time intervals of the 
evolution. The domains of validity of the exponents predicted by perturbation theory, 
-2 near and -3 near z'"*", can be observed here. It is clear that the outermost gridpoints 
in this evolution (using 10000 gridpoints) are indeed located in the "astrophysical zone" 
since r{x - 0.9995) » 2000 3> Amb, where Amb is the Bondi time elapsed since the main 
pulse of radiation has reached the observer at about Ug ~ 3 (see figure 6.8). On the 
other hand, closeness to timelike infinity demands that Awg » r,. Note that r, « r 
for r » 2Mi, where M, « 0.06 is the initial Bondi mass. It is also apparent that for 
Amb ~ ^ the observers are in between the two zones and the power-law exponents seem 
to change smoothly. 

Figure 6.8 displays the power-law decay, with exponent -5, of the Bondi mass for 
the same subcritical evolution. This behavior can be explained by integrating the Bondi 
mass-loss equation (3.7.13) with c oc in the regime of power-law tails. 

Note that in situations involving realistic sources and realistic detectors, power- 
law tails only play a minor (if any) role as an actual signal, but the power-law tail 
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Figure 6.7: This figure shows power-law exponents for a subcritical evolution and illustrates 
the domains of validity of the predictions of perturbation theory for the two zones: -2 near 
and -3 near f . 



results still can provide hints on how to resolve the very important question of whether 
null infinity is a useful idealisation for gravitational wave detectors. Accordingly, we 
suggest to generalize the term "astrophysical zone" to more general, non-perturbative 
situations, by re-enterpreting - in a very loose sense - Awg as a suitably chosen large time 
scale characteristic of the source. The physical idea is that the distance from observers of 
astrophysical phenomena, e.g. gravitational wave detectors, to the radiation sources is 
very large compared to the time during which substantial radiation from the source can 
be observed. This is at least expected for sources where general relativity is important, 
as opposed to problems for which a (Post-) Newtonian approach and the quadrupole 
formula are sufficient. An example would be a binary black hole merger in another 
galaxy, which might have a characteristic dynamical time scale of a fraction of a second, 
and which might be observed during several thousand cycles, including a portion of 
the QNM ringdown. 
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Figure 6.8: We compare the decay of the Bondi mass in supercritical, near-critical and subcritical 
evolutions. In the subcritical case, the Bondi mass is foimd to decay for late times with a power- 
law exponent of approximately -5. 
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Chapter 6: Quasinormal Modes and Tails 



Chapter 




Numerical Evidence for Critical 
Phenomena 

In this chapter we present numerical results of critical behavior in the Einstein-massless 
scalar field system computed using the DICE code and the double-null code. 

Section 7.1 focusses on results from the compactified dice evolution code which 
stress the global observability of critical phenomena at future null infinity. We also show 
results from uncompactified evolutions for cases in which we did not have compactified 
results available. We would like to emphasize that for the local detection of discrete 
self-similarity, inside the past self-similarity horizon, as discussed in section 7.1.3, both 
DICE codes are basically equivalent, the uncompactified DICE code having a slight 
edge over the compactified code as to resolution. 

The results of the double-null code shown in section 7.2 complement the local picture 
given by the DICE code. The double-null code is uncompactified but, in contrast to the 
DICE code, it can penetrate apparent horizons and thus enables us to probe supercritical 
spacetimes a little longer. The reason for introducing this code was mostly to analyze 
the fate of the exterior mass that accumulates outside of the past self-similarity horizon 
in near-critical evolutions, as will be discussed in section 7.1.2. Moreover, it is always 
preferable to have more than one way of computing a numerical result, which, by 
comparison, allows us to rule out numerical artifacts specific to a particular numerical 
scheme. 

7.1 Results from the DICE Code 
7.1.1 Identification of Critical Behavior 

We consider one-parameter families of initial data (p — (pp{uo,x), such that for small 
values of p we have dispersion, while for large values of p we have black hole formation. 
It has been foimd numerically [dT92a, Cho93] that, for any initial data family considered, 
the evolution of near-critical data approach a universal DSS critical solution. In this 
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thesis, we assume universality and restrict ourselves to Gaussian-like initial data 

cpiuo,x)^Ar{xfexp[-('-^^^^fl (7.1.1) 

where r{x) = j^. This choice makes it easy to compare compactified to uncompactified 
evolutions using the same initial data. All results from the compactified DICE code 
presented here use tq = 0.7 and a = 0.3 and a radial resolution of 10000 gridpoints, 
unless stated otherwise. The criticality parameter p is identified with the amplitude A. 

As has become common practice, we find near-critical data through a bisection search 
in p. This procedure yields in particular a numerical approximation to the critical value 
p = p* which defines the threshold of singularity formation. 

In the bisection search, a number of criteria are possible to distinguish dispersion 
from collapse. Their equivalence can be checked a posteriori when comparing the final 
result, i.e. subcritical and supercritical solutions close to criticality. A typical criterion is 
to monitor the ratio 2m /r, where 2m I r - 1 signifies the presence of an apparent horizon. 
This criterion has been used successfully also in combination with slicing conditions, 
that do not penetrate apparent horizons - as is the case in our approach. Remarkably, 
in practice it turns out, that 2m I r > 0.6 is a sufficient criterion to mark a scalar field 
evolution as supercritical. For practical and historical reasons, this is the approach we 
have adopted for our code. 

A number of other options come to mind, in particular in our context of evolving out 
to null infinity, one could e.g. monitor the redshift (2.3.46) or Bondi mass. In a dispersion 
evolution, the redshift will decay to zero, while it will approach infinity when a black 
hole forms. Similarly, the Bondi mass will decay to zero with a characteristic tail 
behavior, as is shown in Fig. 6.8, in a dispersion evolution, and will asymptote to a 
(positive) constant when the field does not disperse. 

The local detection of discrete self-similarity and extraction of the critical echoing 
parameter A ^ 3.44 from a near-critical numerical evolution is discussed in section 7.1.3. 
Here we are concerned with the global picture: In the course of a near-critical evolution, 
remnants of the self-similar dynamics which occur locally, inside the SSH, are radiated 
to future null infinity. Remarkably, we observe that the imprints of DSS behavior are 
still present in asymptotic quantities such as the news function and the Bondi mass (see 
figures 7.1 and 7.8). 

To further test the accuracy of the compactified DICE code in the asymptotic regime, 
we have plotted the left and right hand sides of the Bondi mass-loss equation, (3.7. 12), in 
figure 7.2. Figure 7.3 shows the redshift factor, equation (2.3.46), for the same evolution. 

In addition, in numerical evolutions of supercritical data, the black hole mass has 
been found to exhibit a universal scaling law (5.3.20), in accordance with the discussion 
in section 5.3, (see [HP97, KHA95, Gun97b]): 

In = y ln(p - p*) + W{\n{p - p*)) + const, (7.1.2) 

where y « 0.373 and the function W is periodic with period |A/y in ln(p - p*). Figure 
7.5 shows this scaling behavior of the mass. To analyze the fine-structure of the mass- 
scaling law, one conveniently runs a series of evolutions with masses spaced linearly in 
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Figure 7.1: We show the news function N{u), first as a function of the natural time coordinate 

Mb of an asymptotic observer and also as a function of a suitably adapted time Tb = - In 
where N(tb) is periodic with period A ^ 3.44 after the spacetime has come close to the critical 
solution. Even if the constant m* is not known, it can be determined by a fit to periodicity in Tb. 
Thus, it is possible to observe DSS at J^^ and to extract the critical exponent A. 
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Figure 7.2: This figure superimposes the left and right hand sides (empty and filled circles, 
respectively) of equation (3.7.12), the Bondi mass-loss equation, and plots them against adapted 
time T for a very near critical evolution. The data have been sampled in order that individual 
data points are visible. Until t = 9.0 the empty and filled circles are indistinguishable in the plot, 
i.e. the mass-loss equation holds with a high degree of accuracy. Around x = 11.0 a small black 
hole is starting to form, the resolution has been exhausted and numerical results are unreliable. 
It is apparent that the bursts of scalar field radiation which escape to are periodic in t with 
the same period A/2 which is found in the DSS behavior of fields near the origin. 
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Figure 7.3: This figure shows the redshift factor z = e™ - 1 for the compactified evolution 
depicted in figure 7.2. Variations in the density {2m /r) of the discretely self -similar scalar field 
subject the outgoing light rays to a varying amount of focussing. When a black hole starts to 
form around t = 11.0, the redshift factor diverges exponentially. 



ln(p - p*). One also needs to ensure to measure the mass of the outermost component of 
the horizon in order to obtain the fine-structure of the scaling law shown in figure 7.4. 

This scaling law has been derived by arguments based on linear perturbation theory 
around the critical solution and dimensional analysis however without using a precise 
definition of the black hole mass [ '^QZ, GiinQZl ]. A typical approach in numerical 
simulations based on coordinates that do not penetrate apparent horizons, as far as we 
are aware, is to follow a peak in 2m /r until this quantity almost reaches unity, at which 
point the simulation is usually slowed down by a CFL-type condition, and then read 
off the approximate horizon mass at this point. We follow this heuristic approach, and 
are able to reproduce the mass scaling and fine structure. 

At least conceptually, this approach is problematic without quantifying how much 
mass-energy still remains outside of the peak in 2m /r where the horizon-mass is read off. 
From the perturbation theory argument, one expects scaling behavior only for quantities 
within the self-similar region of spacetime, which does not necessarily extend beyond 
the SSH. We will return to this point in the discussion section 8. 

We have found that the mass scaling law for compactified evolutions exhibits 
a "leveling-off" effect for strongly fine-tuned evolutions. Figure 7.6 shows that for 
In (p - p*) < -30 the logarithm of the detected black hole mass is constant, while further 
away from criticality, the scaling agrees with the expected power law behavior given in 
equation (5.3.20). The fact that this constant is resolution dependent, i.e. it decreases 
with increasing resolution, indicates that it is a numerical artifact of the compactified 
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Figure 7.4: The fine-structure in OTbh after subtracting a linear fit. The measured period 4.6 is 
close to the value predicted by perturbation theory jA/y ^ 4.61. 




ln(p - p*) 



Figure 7.5: This figure shows the mass-scaling law for a series of uncompactified evolutions 
at 5000 gridpoints radial resolution. InniBH is well fit by a straight line with slope y - 0.373. 
For ln(p - p*) < -35 which corresponds to p - p* < 10~^^ our code reaches the precision limit 
for double precision floating point numbers. On the other hand, for ln(p - p*) > -10 we are 
leaving the range of validity of linear perturbation theory which predicts the simple power-law 
behavior of the mass. 
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Figure 7.6: This figure shows the mass-scaling law for a series of 200 compactified evolutions 
using 8000 gridpoints (depicted as a continuous line) and 20 evolutions using 16000 gridpoints 
(depicted as black dots) . Note that as p* depends on the gridresolution N, one has to use p* [8000] 
and p*[16000], respectively for the two datasets. For ln(p - p*) < -30 the black hole masses 
asymptote to a different constant value for each resolution. For ln(p - p*) > -5 we are leaving the 
range of validity of linear perturbation theory which predicts the simple power-law behavior 
of the mass. In the intermediate range -25 < ln(p - p*) < -5, In niBH is well fit by a straight line 
with slope y = 0.373. 



In contrast, uncompactified evolutions generated by the DICE code where the outer 
grid boundary has been adjusted to be close to the ingoing null ray, v = v* , that hits 
the accumulation point, are completely unaffected by this leveling off in the scaling law 
(see figure 7.5). Nor are codes affected that penetrate the apparent horizon and base 
their horizon detection on the first actual apparent horizon detected on a slice, such as 
our double null code (see figure 7.18 in section 7.2.1). 

Note that for a near-critical solution, the approximate value of the accumulation 
time naturally defines an approximate location (i.e. advanced time) of the SSH. 

For observers at null infinity, a natural time coordinate is Bondi time Mb as defined 
in eq. (2.3.35) - Bondi time can be identified with the proper time of timelike observers 
at large distances - see the discussion in Sec. 8. In order to analyze critical phenomena, 
we define a time coordinate which is suitably adapted to self-similar critical collapse 
and set 



The adapted time coordinate Tb can be defined for spacetimes which are close to the 
critical solution inside of the SSH, so that the value of the accumulation time w* can be 
determined by a fit to periodicity in Tb. We have used a fit to periodicity of the news 
function - N(tb) is periodic with period A ^ 3.44 after the spacetime has come close to 



code. 



= -In 




(7.1.3) 
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the critical solution (see figure 7.1). 

Note that Tb is only an approximate adapted coordinate since it depends on the 
relation between Bondi and central time (2.3.35). In order to gain some insight into the 
behavior of Tb(t) consider the simpler case of a CSS collapse. We assume that fi changes 
only little outside the past SSH, i.e. H{u) « f^ssniu). Then, since |Sssh(m) = constant in the 
CSS regime, we find by integrating (2.3.35) that 

Mb ~ Cu, (7.1.4) 

where C = e^^ssH ^j^^j -^^e have chosen the initial condition Wb(0) = 0. Furthermore, it 
follows from the definition of the adapted time coordinates, equations (5.1.7) and (7.1.3), 
that 

Tb « T, (7.1.5) 

in the CSS regime. Numerically, the deviations from this relation for the DSS case turn 
out to be quite small. We cannot expect this relation to hold, if the mass outside the past 
SSH is large compared to the inside, which we observe in the very late stages of strongly 
fine-tuned evolutions. Nevertheless, we have still been able to resolve DSS phenomena 
using the coordinate Tb (see fig. 7.8). 



7.1.2 DSS Behavior in the Bondi Mass and the News Function 

The following argument suggests that the news function is approximately periodic in 
Tb, as shown in figure 7.1. Assume effective DSS data for the scalar field and the metric 
functions on the past self-similarity horizon (SSH) (see figure 7.7). Moreover, we assume 
that the contribution of the right hand side integral in the wave equation for the scalar 
field (3.5.1) can be neglected outside of the SSH, such that the DSS data on the SSH are 
linearly propagated to J^"*" without backscattering, i.e. 

lim ^{u,r) « i/'ssh(w)- (7.1.6) 

Furthermore, assume that changes in |S outside of the SSH are small, so that ^ssn{u) « 
H{u). It then follows that 

n(ub) = 5^ « 

du^ du^ 

= ^[(/)ssh(t)C(t).-m*]|.-2« (7.1.7) 
= e-2feH[^ ((/)ssh(t)C(t)) - (/)ssh(t)C(t)], 

which shows that the news function N{u^) is approximately periodic in t (and in Tb if 
equation (7.1.5) holds) with period A and satisfies N(t + nA/2) = (-1)"N(t). 

In order to determine the behavior of the Bondi mass, we can then rewrite the Bondi 
mass-loss equation (3.7.13) 

dWJn T dUu o 

= -47zN2(tb):^ = -^nule-'^N^iu). (7.1.8) 

«Tb ATb 
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Figure 7.7: A conformal diagram of a critical collapse spacetime. In the backwards light cone 
of the accumulation point the dynamics are close to the DSS critical solution. The lightlike 
boundary is also called the past self -similarity horizon (SSH). Depending on whether the initial 
data are sub- or supercritical, the spacetime will for late times be Minkowski or Schwarzschild. 
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Since N^(tb) is A/2-periodic in Tb, then takes the following form 

"tB(TB) « e-"VssH(TB), (7.1.9) 

where /ssh(tb) is A/2-periodic in Tb. 

This behavior mimics the behavior of the mass-function: We rewrite the mass- 
function in adapted coordinates (t,z), using (5.1.8), 

m{T,z) = ^ze-X{T)u[l - j{z,T)e-^P^^'^^], (7.1.10) 

and evaluate it at the SSH, which, by a judicious choice of the periodic function C,{t) can 
be chosen to be at z = 1 (since the past SSH is a null surface, one needs to ensure that 
VflZ becomes null at z = 1). We obtain 

nissuiT) = e-'f;,^{T), (7.1.11) 

where /s*sh(i^) is periodic with period A/2. 

The periodicity and exponential decay of Wb and mgsH are confirmed by our numerical 
calculations (see fig. 7.8). Note that the Bondi mass levels off at roughly 10~^ of the 
mass in the initial configuration, whereas the mass contained within the backwards light 
cone mssH continues to scale according to the prediction of critical collapse evolution, 
equation (7.1.9). Since the difference, i.e. wJext - m^- ?«ssh/ is almost zero at the initial 
slice, we conjecture that, at late times, wJext is the energy due to backscattering in a critical 
evolution. We do not have a physical explanation for the oscillations seen in ni^xi for 
earlier times. This behavior may, however, be due to errors in the determination of the 
location of the past SSH. In figure 7.9 it is apparent that wJext is localized in a region 
outside the past SSH. 

We observe that ni^xi contained in a slice close to horizon formation is almost constant 
for different near critical evolutions (down to the numerical limit of fine tuning). In the 
very late stages of the evolution, the growing redshift, |S — > co effectively halts the 
numerical evolution, while the error norm max |Euur(TB > 10)| approaches 10~^. 

7.1.3 Local Detection of Discrete Self-Similarity 

We discuss some methods to extract the echoing exponent of discrete self-similarity from 
a near-critical numerical evolution while restricting attention to the past self-similarity 
horizon. 



max(2m/r)-Plots 

The black-hole formation diagnostic 2m I r furnishes a necessary condition for discrete 
self -similarity, namely periodicity in the function 

max — (t), (7.1.12) 
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Figure 7.8: This figure plots the Bondi mass Tfig against both and the adapted time Tb for a 
barely supercritical evolution with final black hole mass Mf ^ 5x 10~^. The Bondi mass ntg and 
the mass at the past SSH, OTssh, are found to decrease exponentially in Tb (with an overlayed Tb- 
periodic oscillation with period A/2), once the evolution has sufficiently approached the critical 
solution near the center of spherical symmetry We also show niExj, the energy present outside 
of the SSH. 
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Figure 7.9: Shown is the mass-function on two null-slices very late in the nearcritical evolution 
depicted in fig. 7.8. Excluding the gridpoint at , the compactified grid extends to large 
values of the areal radial coordinate r. For last slice shown, the location of the past SSH is 
approximately at r ~ 0.0001. The matter exterior to the SSH, is mostly concentrated in the 
radial region 0.5 <r < 100. 
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Figure 7.10: This figure shows a surface plot of |i/'(t,x)| for the same nearcritical evolution as in 
figures 7.8, 6.1 and 6.2. When the initial Gaussian reaches the origin, it is "instantly" (in retarded 
time u) radiated to future null infinity (located at x - 1) by interfering nonlinear ly with 
the field that has not yet reached the origin. Once the evolution has come close to the critical 
solution, the matter field ip{u, x) = (pr decays exponentially; further self -similar features are thus 
not visible in this plot. 

where the maximum is taken over u - const null slices and m{u, r) denotes the mass- 
function in spherical symmetry, as defined in equation (3.7.1). 

The viability of this approach depends on the choice of initial data. For the usual 
Gaussian initial data family, experience shows that the maximum is well-defined. For 
other choices, such as the double-Gaussian initial data of section 7.2.4, two (or more) 
peaks in 2m jr may compete with each other and thus the maximum can jump from one 
slice to another. 

The echoing exponent A can be determined by analyzing plots of max(2m/r) vs. t 
and adjusting u* such that successive echoes obey a constant spacing, which equals A/2. 
In practice, a crude estimate of u* can be obtained in slightly supercritical evolutions 
as the coordinate time of the formation of the (small) black hole. This process can be 
carried out by hand or automatized by a fitting procedure. The initial transient in which 
the system first approaches the vicinity of the attractor has to be ignored. The result of 
such a fit is shown in figure 7.11. 

DSS-Snapshots 

Although plots of the diagnostic max(2m/r) are handy for determining A and u* , since 
they are easy to generate, they do not encompass the full scope of DSS, since, by taking 
the maximum, the information available is collapsed down to just one data point per 
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Figure 7.11: This is a typical plot of the black hole diagnostic max2m/r vs. t computed using 
the compactified DICE code with 8000 gridpoints. The parameter A can be readily extracted 
from the picture. The accumulation time u* was determined as the time of the formation of the 
very small black hole with mass wZbh - 3.5 X 10~^. For t > 12, max2?n/r rapidly approaches 
unity, which heralds the formation of an apparent horizon. 



null-slice. Given a function on spacetime which is DSS, it shows this behavior at 
any given radial position. This much stronger criterion can be investigated with the 
following method (adapted from [Cho98]). 

In coordinates {t, r) DSS corresponds to a logarithmic rescaling in both time and 
space. As has been discussed in section 5.1, in adapted coordinates (t,z) DSS amounts 
to 

(/)(T + nA/2,2) = (-ir(/)(T,z). 



Here, we prefer to use In r instead of z. From the definition of z. 



z - 



C(T)t/»' 



observe that 



Inz = Inr + T - ln(C(T)M*). 
Thus, the following equality holds 



(7.1.13) 



(7.1.14) 



(7.1.15) 



(/)(t + n A/2, Inr + t + nA/2 - ln(C(T + nA/2)u)) 

= (-1)" (/)(t. In r + t - ln(C(T)M*)), (7.1.16) 
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and, for a given t = const, and using C,{t + nA/2) = Q{t), we have 

(/)(T + nA/2,lnr + nA/2) = (-1)" (/)(T,lnr). (7.1.17) 

Successive frames in the plots shown in figures 7.12 and 7.13 are equally spaced in 
the adapted time coordinate t. The quantities are plotted against In r with appropriate 
"shifts" in the null slices according to equation (7.1.17). 

While dimensionless quantities {(p, 2m/r, ^, y) are mapped according to (7.1.17) 
without a change of scale under the diffeomorphism associated with DSS, dimensionful 
quantities such as or also scale according to their dimensions, e.g. [^] = L and thus 

i/^(t + A/2, In r + A/2) = -e^/^i/;(T,lnr). (7.1.18) 

As an example, figure 7.14 shows the scalar field and its image under the DSS diffeo- 
morphism. 

To be able to more clearly visualize discrete self-similarity, we show the scalar field 
as a function of t and - In x for a near-critical compactified evolution in figures 7.15 and 
7.16. Figure 7.17 depicts the ingoing null geodesic that hits the accumulation point of 
the DSS-solution. 



Results from the DICE Code 



117 



T = 2.16 


T = 2.60 


T = 3.03 


T = 3.47 










X = 3.90 




T 4.77 


T = 5.20 


— — ""^^^ 






x = 5.63 


T = 6.07 


T 6.50 


T = 6.94 








T = 7.36 


^^7^0 


T = 8.23 


T = 8.67 











Figure 7.12: This figure displays snapshots of a near-critical uncompactified evolution with 
5000 gridpoints of the massless scalar field (p as a function of ]n r. The frames are evenly spaced 
in the adapted time coordinate t = - In ^^-^ . Observe that the shape of cp is identical in every 
8th frame; thus, cp is periodic with period A ^ 3.44. Moreover, (p is antisymmetric with respect 
to every 4th frame, i.e. it satisfies the half-period self-similarity condition (5.1.23). Figures 7.13 
and 7.14 show other discretely self -similar quantities from the same numerical evolution. 



118 I Chapter 7: Numerical Evidence for Critical Phenomena 



T = 2.16 


T = 2.60 ^ 


^^^T=^^^03^^ 


X = 3.47 




T = 4.33 


X = 4.77 ^ 


X 5.20 


T = 5.63 


X = 6.07 


X = 6.50 ^ 


X = 6.94 


T - 7.36 


T = 7.80 


x = 8.23 ^ 


X = 8.67 



Figure 7.13: Snapshots of the dimensionless black hole formation diagnostic 2m/r for the same 
near-critical evolution as shown in figure 7.12. 
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Figure 7.14: This figure shows i/'(To,lnr) (depicted by a few sampled dots) at a time tq = 
4.77, when the evolution was in the echoing region overlaid with its image under the DSS 
diffeomorphism Oa, e^ip{TQ + A, Lnr + A) (represented by the continuous line). Since is not 
dimensionless, its amplitude had to be rescaled by e^. 
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Figure 7.15: This figure shows |(^(t, -lnx)| for a near- critical compactified evolution. In the 
direction of increasing t the features of the scalar field repeat themselves with a period A/2 
since (p fulfills the half-period self-similarity condition (5.1.23). Concurrently, these features are 
shifted towards the origin by A/2 in - In r (In x ~]nr near the origin) . Figure 7.16 shows a planar 
view of the same surface. 
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Figure 7.16: This figure shows a planar view of |c/)(t, - lnx)| for the same evolution depicted 
in figure 7.15. The ingoing wave packet is clearly visible on the initial slice (t - 0) of the 
evolution. The the tail ends of the bumps are enveloped by the null geodesic that barely hits the 
accumulation point which is shown separately in figure 7.17. 
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Figure 7.17: This figure shows the null geodesic which almost hits the accumulation point in 
the evolution shown in figures 7.16 and 7.15 
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7.2 Results from the Double-Null Code 



The results shown here are based on fine-tuning of Gaussian initial data 

(7.2.1) 

In section 7.2.4 we introduce double Gaussian data that are merely a superposition of two 
such Gaussians with parameters adjustable for each pulse. 

7.2.1 Scaling 

To complement the results from the DICE code, figures 7.18 and 7.19 show the ex- 
pected scaling behavior for the black hole mass and the Ricci scalar curvature from a 
one-parameter family of evolutions generated by the double null code, respectively. 
According to its dimension being 1/L^, the scalar curvature scales as (p* - p)~^'', as 
discussed in section 5.3. 



s{u = 0,v) = A exp 




-40 -35 -30 -25 -20 -15 -10 



ln{p — p*) 

Figure 7.18: This figure shows mass scaling for a family of evolutions using 8000 gridpoints. 
The critical exponent is given by the slope of the fitted line, which yields y ~ 0.373. Note that 
this scaling is perfectly consistent with the analytical scaling law and does not suffer from the 
leveling-off artifact present in mass scaling laws generated from the compactified DICE code. 



124 I Chapter 7: Numerical Evidence for Critical Phenomena 



I 



35 



30 



25 



20 



15 - 



10 




-40 -35 -30 -25 -20 -15 -10 -5 



ln|p-p*| 

Figure 7.19: Here we show curvature scaling for a series of subcritical evolutions using 8000 
gridpoints. The maximum of the Ricci scalar at the axis, plotted as b-i(l - R), behaves as in 
equation (5.3.12). The scaling exponent extracted from the fit is y ~ 0.374. 



7.2.2 Horizons and m^xr 

In section 7.1.2 we mentioned that exterior matter, rngxT/ accumulates outside the past 
self-similarity horizon and we conjectured that, for late times, the mass originates from 
backscattering of outgoing radiation during the self-similar collapse. In a supercritical 
evolution, two scenarios are possible: Either niEXT eventually falls through the horizon, 
or it is radiated to infinity. In the former case, the resulting black hole will have a tiny 
but finite Bondi mass, no matter how fine tuned the data are, whereas, in the latter case, 
the Bondi mass can be made arbitrarily small. 

We would like to explicitly verify whether this matter, tn^xj, falls through the event 
horizon. (A compactified code that could penetrate apparent horizons would presum- 
ably still have to stop at the event horizon, due to reaching future timelike infinity, 
Using an uncompactified code, we will have to settle for a reasonable local approxima- 
tion of the event horizon.) 

At first, it would seem that one could investigate the final fate of tn^xj with any 
evolution code that can penetrate apparent horizons, such as the double null code of 
section 4.2. As it turns out, however, codes based on outgoing null slices are ill-suited 
to tackle this issue. 

Only the immediate vicinity of the event horizon is relevant for the solution of 
this problem. Null slices too far to the future of the event horizon no longer contain 



Results from the Double-Null Code 



125 



wJext/ since the slices will have "bent" towards the curvature singularity at r = 0, and 
trapped gridpoints which come too close to the singularity will have been excised, as 
has been discussed in section 4.2.4. This leaves us with a slice where the very features 
we would like to study, namely the exterior matter, have probably already fallen into 
the singularity. 

To illustrate the behavior of the mass function on outgoing null slices in the vicinity 
of apparent horizon formation we compare a series of slices for a near-critical evolution 
using 5000 gridpoints in figure 7.20. If a slice contains an apparent horizon, the function 
r{v, u = const) increases monotonically from the origin outwards until the locus of the 
AH is reached; there, it reverses direction and decreases monotonically until it reaches 
the singularity at r = 0. Gridpoints which come too close to the singularity must be 
excised. In contrast, figure 7.21 depicts a series of slices for a near-critical evolution 
which ultimately is subcritical. The two plots are superimposed in figure 7.22. 

To improve the chance to observe the exterior mass falling into the horizon, it is 
clearly beneficial to decrease the timestep Au as far as possible. Unfortunately, the 
double null code relies on the condition Au - Av for the approximation of derivatives 
close to the origin. We have tried to use interpolation to get rid of this condition, but our 
efforts introduced instabilities that made evolutions unusable. There is, in principle, 
another way. Namely, keep the condition Au — Av and increase the overall spacetime 
resolution, i.e. both gridspacings. The result of such a brute force evolution using 
100000 gridpoints is shown in figure 7.23. The horizon mass is almost two orders of 
magnitude lower, but the bending-over of the slices is still much too abrupt and all the 
exterior mass has been excised on the first slice which contains an apparent horizon in 
this evolution. 

The behaviour discussed above is a deficiency intrinsic to our choice of slicing: 
Null-slices are inherently unstable near a horizon. Outgoing null-slices either extend 
till J^""", if they are untrapped, or end at the singularity r = (with the exception of the 
event horizon which extends till /"*"). Thus, it would be seem to be appropriate to use a 
Cauchy code to investigate this issue further. This would allow us to scan the horizon 
for exterior infalling matter. 



7.2.3 Dependency of niext on Initial Data 

In the following, we would like to discuss the dependance of the exterior matter, m^yj, 
on the initial data. First, choose a family of initial data and fix all free parameters except 
for the parameter p that will be used for fine-tuning. Within such a family, as we tune 
the parameter p to criticality we can reliably measure Wext only if at some null slice 
in a given evolution, the interior mass is of the same order or smaller than /Wext or the 
interior and exterior mass are spatially well separated. In practise, this condition is 
only fulfilled at late times for very near-critical evolutions which severely restricts the 
interval in the initial data parameter and evolution time over which we can accurately 
measure the exterior mass mExi- 

From the analysis of a typical single near-critical (but supercritical) evolution, m^xi 
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Figure 7.20: Shown is the mass-function of a series of null slices of a supercritical evolution 
close to criticality. The formation of an apparent horizon manifests itself in the bending over of 
the curves in the r-coordinate. After having reached its maximum r-value the slice bends back 
towards r = and reaches the curvature singularity. This necessitates the excision of gridpoints 
coming too close to r = 0. 
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Figure 7.21: Shown is the mass-function of a series of slices of a subcritical evolution close 
to criticality. The initial data are very close to those in figure 7.20, which ultimately forms an 
apparent horizon. In this evolution, the scalar field disperses, ultimately, reaching Minkowski 
space 
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Figure 7.22: This figure shows figures 7.20 and 7.21 superimposed. The two evolutions almost 
coincide in the first few slices shown, when the apparent horizon has not yet formed. Gradually, 
the agreement gets worse for r < 10~*, while it is still close for larger r. Then, an apparent 
horizon forms in the supercritical evolution, while the subcritical evolution continues to radiate 
scalar field that has reached the origin and ultimately disperses. 
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Figure 7.23: This figure shows the mass-function before and after the formation of an apparent 
horizon for a brute force evolution at 100000 gridpoints. The reduction of Am was still not able 
to resolve the bending over of the slices fast enough to have the trapped slice extend up to the 
exterior matter. 
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oscillates from the time we can first detect it until the last null slice before horizon 
formation, while the interior mass (contained in the past SSH) undergoes critical collapse 
and periodically sheds mass. If we instead consider a family of supercritical evolutions 
which become more and more near-critical, and measure m^xr at the last slice before 
horizon formation, m^xi roughly stays constant. 

It is also important to note that, using generic (e.g. Gaussian) initial data, which have 
been fine-tuned, the solution evolves through a non-universal and family dependent 
stage before it approaches the self-similar regime. Since backscattering of radiation is 
already present in this initial stage, we expect mExr to also depend on the shape or family 
of initial data used. 

7.2.4 Double Gaussian Initial Data Results 

In an effort to shed some light on the issue concerning the possible infall of the external 
mass niEXj, we have devised initial data to model this physical situation. Essentially, we 
are dealing with a near-critical evolution where the interior part is in the DSS regime (for 
a certain time) while further outside, there is another concentration of energy. Therefore, 
it is sensible to chose initial data which consists of two Gaussian pulses of scalar field, 
which we call double Gaussian initial data for short. 

As it turns out, Wext is not dense enough to form an apparent horizon by itself. This 
can be guessed from the fact that there is no apparent peak in 2m/r at the locus of Wext- 
It is also possible to evolve data that are based on the numerical scalar field solution 
on a null slice of a near-critical evolution where the presence of Wext is the dominating 
contribution of mass. Using such a solution, s(u, r) = s{u,v{u,r)), (preferably the last 
available slice before the formation of an AH) as new initial data so(m = 0,v) = s{u, 2r) is 
fairly straightforward. In is necessary to employ interpolation to obtain equally spaced 
initial data so(0, v) for the scalar field . Here, we have also had to fix the null gauge, 
which determines how the gridpoints are labeled by u and v. If necessary, taper the 
gridfimction such that the gradients in the interior part are smooth enough that they 
do not form trapped surfaces; in practise, this is taken care of automatically by the 
interpolation. We have carried out this procedure and found that the new spacetime is 
subcritical. 

The accumulation of exterior mass that occurs naturally via backscattering in critical 

collapse could be modeled by data where the outer pulse is chosen weak enough so 
that it does not collapse by itself. Unfortunately, simulations of that type did not yield 
any new insights, since we still run into the problem of the bending-over of slices and 
forced excision of the part of the slice we would like to follow further. 

Although not directly related to the issue we would like to solve, we have inves- 
tigated spacetimes where the outer pulse can form a horizon, while the interior pulse 
evolves, more or less unaffected from the outer pulse, as they are interesting in their 
own right. In particular, we consider a configuration where the outer pulse froms an 
AH first (in retarded time u), while we can still follow the interior pulse in its critical 
collapse. Later (in retarded time), the interior pulse may itself form trapped surfaces. 



130 



Chapter 7: Numerical Evidence for Critical Phenomena 



which creates a jump in the trapping horizon. A spacetime diagram of such an evolution 
is shown in figure 7.24. Figure 7.25 presents a density plot of 2m/r for this evolution, 
(with a detail shown in figure 7.26), while figure 7.27 shows a perspective view. 

7,2,5 Timelike Observers and t versus r Diagrams 

For timelike observers at r = const we can compute the proper time t{u, r - const) as 
discussed in sections 3.8.3 and 4.2.7. In doing so, we obtain a {t{u, r), r) grid on which 
we can display quantities of interest, such as the density function 2m jr. Since this grid 
is neither orthogonal, nor equispaced, we need to resort to interpolation for some plots. 

In order to have more resolution in the vicinity of the origin, one can use a static 
mapping in the distribution of the r — const observers. A useful convex mapping is to 
equally space the observers in a computational coordinate E, and define r{E,) = L{e''^ — 1) 

with r{E, = fmax) = ''max/ SO that L = rmax/(e'"''"^'' - 1). Figure 7.28 shows a color-coded 
density plot of 2m I r on such a grid, where the null slices u = const have been overlaid. 
Grid refinements are clearly visible when the density of the slices suddenly doubles. 
The geometry is initially flat. The evolution uses double Gaussian initial data, where 
the exterior pulse forms the first (in retarded time u) apparent horizon and the interior 
collapses at a later null slice after moving away from the DSS regime. 

In contrast, figure 7.29 compares the fate of r = const observers in a supercritical and 
a subcritical evolution. 
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Figure 7.24: Similar to figure 4.9, this spacetime diagram shows the formation of an apparent 
horizon and the region of spacetime that has been evolved. This time, we are dealing with 
an evolution using double Gaussian initial data, chosen such that the outer pulse forms the 
apparent horizon (first in u) while the inner one reaches 2m /r - 1 later in retarded time u. Then, 
excision kicks in and we are left with a very small grid. 
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Figure 7.27: This figure shows a perspective view of the spacetime depicted in figure 7.25. 
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Figure 7.28: This figure shows a density plot of the diagnostic 2m I r on a {t,r) grid for a 
supercritical evolution using double Gaussian initial data. 
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Figure 7.29: This figure shows spacetime diagrams for r = const observers. In the left evolution, 
the scalar field configuration formed an apparent horizon as is evident from the r = const curves 
approaching the null surface r = 2M of the event horizon. In contrast, the evolution shown 
in the right frame turned out to be subcritical. Compare this figure with the Penrose diagram 
shown in 1.1 
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Chapter 

Discussion 



In this thesis we have presented numerical constructions of portions of near-critical 
spherically symmetric spacetimes that extend up to future null infinity. The simulations 
are based on a compactified approach, where the equations have been regularized 
in a neighborhood of null infinity by introducing the mass-function as an additional 
independent variable. Further simulations have been carried out using a double-null 
code that can penetrate apparent horizons. In this discussion, however, we focus on key 
results obtained from the compactified code as they represent the new developments 
in this thesis. The resolution necessary to resolve critical phenomena is gained through 
letting our gridpoints fall along ingoing null geodesies, following Garfinkle [ ] . The 
grid is repopulated with grid points, which are filled with values through interpolations 
when half of the gridpoints have reached the origin. This is not optimal, but ensures 
sufficient resolution for our current purposes. 

We reproduce standard features of near-critical solutions, such as "echoing" and 
mass scaling with fine structure. The mass scaling law for compactified evolutions 
shows a "leveling-off" effect for very near-critical solutions that is not present for un- 
compactified simulations which we attribute to numerical errors. In addition, we extract 
radiation at null infinity by computing the news function and find a signal with rapidly 
increasing frequency as measured in Bondi time, which is the natural time coordinate 
associated with far away observers. In order to simplify the analysis of near-critical 
spacetimes, we have defined an approximate adapted time in terms of Bondi time, in 
analogy to the standard time coordinate which is adapted to the DSS discrete diffeomor- 
phism (5.1.1). This coordinate can be used by asymptotic observers to render the signal 
from a near-critical collapse (almost) periodic. The fact that such a definition actually 
works out, and makes DSS periodicity manifest in quantities defined at null infinity is 
a non-trivial result, which we explain in Sec. 7.1.1. 

Note also that the amplitude of the news function, shown in Fig. 7.1 stays fairly 
constant after the initial transient. This feature is clearly universal, as long as the 
radiation signal is dominated by the DSS collapse, since the system can be approximated 
in the DSS region by a perturbation of the critical solution. Thus, in section 7.1.2, we 
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neglect contributions from scalar field which is far outside the DSS region and does not 
contribute to the critical collapse dynamics. Consequently, in this scenario the essential 
free parameters determining the radiation signal from an actual near critical solution 
are the number of cycles the solution spends in the neighborhood of the critical solution 
and the length scale (e.g. the "size" of the past SSH) at which the solution comes close 
to the critical solution. The robustness of this scenario, i.e. what happens if the initial 
data are such that there is significant mass outside of the SSH, is beyond the scope of 
this thesis and an issue for future research. 

Perhaps the most surprising feature of the radiation signal has emerged from our 
investigation of QNMs, which has been motivated by [GPP941)], where the first quasi- 
normal mode is found in collapse evolutions and the question is posed as to how QN 
ringing would change close to criticality. We find that even in very close-to-critical 
evolutions there is a correlation of the radiation signal with the period of the first 
quasinormal mode, determined from the time-dependent value of the Bondi mass, as 
discussed in Sec. 6.2, Figs. 6.1 and 6.2. This correlation between the radiation of the 
highly dynamical near-critical solution and the quasinormal mode, that is defined in 
terms of perturbations of a static spacetime, certainly deserves further investigation. 
This surprising feature might even turn out to be a key toward understanding the phe- 
nomenon of DSS behaviour in near-critical spacetimes. Our results seem to suggest 
that the effective curvature potential for a DSS self-gravitating field acts as a quasi- 
stationary background for scattering processes which can be approximately described 
by quasinormal modes of a one-parameter family of Schwarzschild black holes with 
exponentially decreasing mass. 

This Schwarzschild background with decreasing mass can also be modeled by the 
Vaidya solution as discussed in section 6.3. Abdalla, Chirenti, and Saa [ACS06] have 
computed QNM for the Vaidya metric and have put forth a criterion to analyze whether 
a spacetime is in the "stationary adiabatic regime" where the real part of the QNM 
varies inversely with the mass function. We have been able to verify that our critical 
collapse evolutions are indeed in this stationary regime. 

The question of the applicability of QNM-motivated estimates is quite relevant for 
numerical relativity, e.g. when extracting wave forms from binary mergers. 

For subcritical initial data we can evolve for very long times, and thus are able to 
observe power-law tail behavior as shown in Figs. 6.7 and 6.8. Analytical calculations 
predict different falloff rates for radiation along null infinity and along timelike lines, 
and the natural question arises, which falloff rate would be seen by a hypothetical 
observer (in a realistic case, observation of power-law tails would require an extremely 
large signal-to-noise ratio). Accordingly, our results depicted in Fig. 6.7 show how the 
rates at finite but large radius correspond to the value for null infinity for a while before 
they approach the expected late-time value for finite radius. The interpretation of this 
phenomenon is suggested by perturbative work of [BarQQa, Bar99b], where different 
tail falloff rates are computed. There, different regions of spacetime are identified, 
where certain approximations hold. Within the perturbative regime, results obtained 
for null infinity are valid for what has been termed the "astrophysical zone" by Leaver 
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[Lea86a, Lea86b, Bar99a], and which is defined as the region where Aug r. (Note, 
however, that Ub » M must also be satisfied). The physical idea is that the distance 
from observers of astrophysical phenomena, e.g. gravitational wave detectors, to the 
radiation sources is very large compared to the time during which the source radiates 
at an observable rate. 

We argue that our results illustrate that the relevant falloff from the point of view 
of an astrophysical observer is the falloff rate at null infinity, in accordance with the 
prediction from perturbation theory. We believe that this is a nice model calculation 
that exemplifies that null infinity is indeed a useful approximation for observers at 
large distance from the source, in the sense that such observers are located in the 
"astrophysical zone". 

Along the same lines, we would like to point out that by taking appropriate limits 
in a conformally compactified manifold, worldlines of increasingly distant geodesic 
observers converge to null geodesic generators of future null infinity and proper time 
converges to Bondi time [FraOO]. Note also that a naive correspondence between ob- 
servers at large distance and spatial infinity would be problematic, e.g. compactification 
at spatial infinity leads to "piling up" of waves, whereas at null infinity this effect does 
not appear - waves leave the physical spacetime through the boundary at null infinity. 

Under practical circumstances, null infinity more realistically corresponds to an 
observer that is sufficiently far away from the source to treat the radiation linearly, but 
not so far away that cosmological effects have to be taken into account. We are however 
not aware of a discussion where this sloppy picture has been made more precise. 

When looking at critical collapse from a global spacetime perspective, as we have 
done here, one is confronted with some issues concerning the mass-scaling, that we 
would like to comment on briefly: When asking the question of whether infinitesimally 
small black holes can be formed - the question which triggered the original work on 
critical collapse - it could be phrased in two slightly different ways: (i) can we form 
black-hole solutions where the final-state black hole has arbitrarily small mass, or (ii) can 
we form arbitrarily small apparent horizons. The question which has been answered in 
the affirmative by critical collapse research is the second one. The first one still seems 
open. Our results for near-supercritical evolutions indicate that the mass outside the 
past self -similarity horizon, is significantly larger than the mass leading to scaling. We 
conjecture that this external mass originates from backscattering of outgoing radiation. 
If this mass falls into the black hole, then, clearly, one cannot form arbitrarily small black 
holes, no matter how fine-tuned the data are. 

The application of a double-null code that can penetrate apparent horizons has not 
been able to decide whether the exterior matter really falls through the horizon; instead 
of a null slicing, which has proven inherently unstable close to the event horizon, we 
propose to use a Cauchy evolution code for future work. 

Considering the observability of critical collapse from future null infinity in an 
astrophysical context, one might first object that the matter model discussed here is 
unphysical. On the one hand, it was chosen for mathematical simplicity. On the other 
hand, the curved space wave equation of the simple scalar field model still shares some 
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essential nonlinear features inherent in the full Einstein equations without symmetry. 
Moreover, it would be possible to remedy this deficiency by using e.g. a perfect fluid or 
gravitational waves in axisymmetry. A much more serious problem is the high degree 
of finetunedness required in the initial data to develop sufficient self-similar features in 
spacetime such that they could be detected by far away observer. 
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Tensor Components in Bondi 
Coordinates 
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Energy momentum tensor of a massless scalar field 

Tab = V«(^)Vi,(^) - ^g„i,Vc(^)V^(^) (A.0.14) 

Tun = V + ^7 (7(^')' - 2^^') (A-0-15) 

Tur = \^W? (A.0.16) 

Trr =Wf (A.0.17) 
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T^^ = sin2 \^{<^'f - 24)(^)') (A.0.19) 

(A.0.20) 

T = g^'Tab = -e-^^ ^-Wf - 2#') = ^«T^B (A.0.21) 
Curved space wave equation (matter field equation) 

= = [(^ + (-)')c|)' - -c/) - 2c|)' + -c|,"J (A.0.22) 
Christoffel symbols 

= ^ - ^''^iS' - rV + 4r2|S) (A.0.23) 
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Tensor Components in Double Null 
Coordinates 



Einstein Tensor 
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Christoffel symbols 
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Appendix 

Numerical Methods 



C.l Finite Difference Methods for ODEs 

Following [AP98], we consider numerical solutions to the initial value problem (IVP) 

dy 

^ f{t,y{t)), yito) = yo, 0<t<T (C.1.1) 

on a mesh 

= < ^1 < • • • < tN-i <tN^T (C.1.2) 

andlet hfi = tn—t„-^ the nth step size. We assume sufficient smoothness and boundedness 
on f{t, y) so that a unique solution y{t) exists and has as many bounded derivatives as 
required. Let i/;, be the mesh function which takes on the value i/„ at each f,„ where the 
yn denote approximations to the exact solution y{tn). 

For the sake of simplicity, consider the forward Euler method, 

y^^' = y^ + hn^if{t\y^). (C.1.3) 

We define the difference operator Nh of the forward Euler method as 

Nnu{tn) := ~ "(^"-i) _ f^tn-l, U{tn-l)), (C.1.4) 

for some function u{t) defined at the mesh points. Then, Njiyh = represents the 
discretized form of the ODE (C.1.1). The local truncation error (LTE) (or local discretization 
error) ^ is the residual of the difference operator when it is applied to the exact solution: 

Tn := Nhy{t„). (C.1.5) 

' This definition of the local truncation error can also be interpreted as the error generated by a single 
numerical step starting from exact data and dividing the result by h,,: i.e. for the forward Euler method 
denote the result of this step by y„ := y{t„-i) + /2„/(J„^i,i/(t„-i)). Then t„ = {y{t„) - y„) /h„. In contrast to 
this definition, the LTE is often defined by i/(f„) - y„. 
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The difference method is said to be consistent (or accurate) of order p if 



Tn = 0{K) (C.1.6) 

for a positive integer p. For the forward Euler method we have 

Tn = ^^^"^ ~ J^^""'^ - f{tn-l, y{tn-l)) = y + OQll), (C.1.7) 

thus, it is consistent of order 1. 
We define the global error 

£>„ := yn - y{t„), (C.1.8) 

and assume that eo = 0- Let H - maxi<„<]v h„ and assume that NH is bounded indepen- 
dent of N. The difference method is said to be convergent of order p if the global error, e„, 
satisfies 

en = Om, (C.1.9) 

forn = 1,2,..., N. 

Theorem 3.1 of [AP98] states that: If a numerical method is consistent of order p and 
fulfills a certain concept of stability ("0-stability" ^ ), then it is convergent of order p. 

In contrast to the local truncation error, the local error is defined as the amount by 
which the numerical solution i/„ at each timestep differs from the solution y{tn) to the 
initial value problem 

yit)=fit,yit)) yitn-i) = yn-1. (C.1.10) 

The local error is defined by 

In = y{tn) - Vn- (C.1.11) 

For the numerical methods we consider here, the local error indicators h„T„ and /„ are 
closely related: 

hn (|t„| + 0(HP+1)) = U (1 + 0{hn)) . (C.1.12) 



C.2 Runge Kutta Methods 

We discuss the Runge-Kutta methods [AP98] used to solve ODEs in our numerical codes. 
Runge-Kutta methods belong to the class of one-step methods, i.e methods that do not 
use any information from previous steps. In a typical step of size h = tn+i - tn, we seek 
an approximation i/„+i to y{t„+i), given the result of the previous step, y„. 

^ A difference method is called 0-stnble if there exist positive constants hg and K such that for any mesh 
fimctions X;, and Z;, with h < hg, 

\<K\\xo-Zo\ + maxlNiMtj) - Mz;,(t,)| }, l<n<N. 

1 l<i<N I 
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Given the scalar initial value problem 

^=f{t,y{t)), y{tQ)^yo, 0<t<T (C.2.1) 

we can approximate the area under the curve y(t) by applying one of the quadrature 
rules given in section C.2.1 to the integral in 

yitn+i) - yitn) = f y{t)dt (C.2.2) 

on a mesh 

Q = tQ<ti<---< tN-i <tN^T (C.2.3) 

and let h„ = t„- tn-i the nth step size. For instance, we can use the height at the midpoint 
of the interval, i.e. y{tn+iii), where tn+iji = in + hjl. This choice leads to the implicit 
midpoint method 

yn+l =y„+hf [tn+l/2, ^" Y"^' ) • (^-^-^^ 

Approximating y{tn+iii) by the forward Euler method we obtain the explicit midpoint 
method: 

h 

Vn+m = Vn + l^f{tn, Vn) (C.2.5) 
Vn+l =yn+hf {tn+1/2, Vn+lll) ■ (C.2.6) 

According to equation (C.1.6), both midpoint methods have a local truncation error of 
OQi^) and are consistent of order 2. 

Applying the trapezoidal rule to C.2.2 yields the implicit trapezoidal method 

h 

+ ^ [fiin, Vn) + fitn+1, Vn+l)] ■ (C.2.7) 

To obtain an explicit method, we approximate using the forward Euler method, 
which yields 

y„+i = yn + hf{t„, yn) (C.2.8) 
h 

+ o [/(^"' Vn) + f{in+l, Vn+l)] • (C.2.9) 

This is an explicit 2-stage method of order 2, known as the explicit trapezoidal method. 
C.2.1 Quadrature Rules 

For the sake of reference we list some very well known approximations to the integral 

rh 

I = Jg f{x)dx and their global errors: 
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Midpoint-Rectangle rule: 



I = hY^fi.l/2+0{h^) 



i=l 



Trapezoidal rule: 



I = h 



1 1 
-/o+/i + ---+/„_i + -/„ 



1/3 Simpson's rule: 

For the subinterval [i - 1, z + 1]: 



while for the entire interval 



(C.2.10) 



+ 0{h^) (C.2.n) 



J^' f{x)dx = ^ + 4/, + fi^t] + 6)(;25), (C.2.12) 



^ = 3 [/o + 4/i + 2/2 + • • • + 4/^-1 + /,„] + ^(/z^)^ 



(C.2.13) 



where m is even. 



C.3 Error Analysis of the NSWE-algorithm 
C.3.1 Local Errors 

The local errors in the NSWE-scheme (3.5.2) of section 3.5 stem from the approximation 
of the integral by evaluating the integrand at the center of the small null parallelogram 
(see figure C.l) and from the multiplication with the area (Au)(Ar) of this diamond- 
shaped region. Evaluating the integrand at the center amounts to the same as using the 
average between the points E and W. 

Assume that we know a function / at the endpoints i]o and r]i of a given interval 
(see figure C.2) and want to evaluate it at the midpoint of the interval r]mid- Using linear 
interpolation at the center of an interval [rjo, rji] we find 

/mid «^[/(?]0) +/(?]!)]. (C.3.1) 

The approximation error for the linear interpolation polynomial with pairwise distinct 
nodes is (see [DH03]) 

/(^mid) - P(/h0, ^]l)(?]mid) = n " " (C.3.2) 

for some £ irjo,rji). Therefore, this polynomial approximation is second order accurate. 



Error Analysis of the NSWE-algorithm | 147 



N 




S 



Figure C.l: The diamond-shaped region in the NSWE-scheme. 

h h 

r V ^ ^ 

^0 rjmid ^1 

Figure C.2: Illustration for the midpoint rule. 

The approximation of the integral in (3.5.2) applies the midpoint-rule to the integrand 
which is multiplied by the area {Au){Ar) which is of the order 0{h^). (If Au « Ar, then 
(AM)(Ar) « 2h^.) Since both factors are second order accurate, the local error of one 
NSWE-integration step is 0(h'^), assuming that il>w, 4'e and ips are known at least to this 
order. 

C.3.2 Global Errors 

In the following, we give a heuristic argument for the global error of the NSWE-scheme. 
First, we restrict analysis to the fixed background case. This entails that the NSWE- 
scheme is no longer coupled to the geometry and null-geodesic ODEs, which would in 
general have to be solved concurrently. To simplify matters even more, we discard the 
Taylor expansion-region and use just ip(r = 0) = as a startup condition. Note that for 

flat space the algorithm is exact, since (y) =0. 

The general idea is the following: We start out from the initial slice where ip is 
known exactly. On the next slice, we march outwards from the origin using the NSWE- 
algorithm whose computational molecule is shown in figure C.3 and accumulate local 
error terms of order 0{h'^) as illustrated in figure C.4. If we make n steps and /i « 1/n we 
will loose one order of accuracy. On the next slice, however, we are still using i/^-values of 
order 0{h^) near the origin. In a Cauchy code we would accumulate errors only along 
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N 




Figure C.3: The computational molecule for the NSWE-scheme. 



the time direction. Here, we can distinguish two directions of error accumulation: 
one radially outward and another along the ingoing null-geodesics. Thus, we loose 
two orders in total and the diamond algorithm turns out to be globally second order 
accurate: 

4> = 4>true+Oih^). (C.3.3) 

For the coupled NSWE-algorithm one has to choose the orders of accuracy for the 
ODE solvers such that the diamond scheme stays globally second order accurate. In 
practise, second order accurate schemes (see C.2) have proven sufficient to ensure this. 
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NSWE-scheme Cauchy-scheme 

Figure C.5: This figure shows the directions of error accumulation present in the NSWE-scheme 
and a protot5^ical Cauchy finite differencing scheme. In the Cauchy scheme, we loose only one 
order of accuracy from the accumulation in the time direction. In the NSWE scheme, however, 
there is coupling in the radial direction, as we see from its computational molecule (figure C.3). 
Combined with the ingoing null direction we therefore loose two orders of accuracy, in total. 
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D 



Convergence Test Methodology 



We briefly describe how we can check the accuracy of a numerical evolution code using 
the concept of convergence tests. These tests obtain error estimates by subtracting 
numerical solutions on a grid which is shared by a set of evolutions using different 
(usually by powers of two) gridspacings. Depending on whether the quantities to be 
tested are zero on the continuum level or not, we take one of the following approaches: 
In a so-called 3-level convergence test, we use three different grid resolutions, with n 
being the number of gridpoints in the lowest resolution grid, and the higher ones having 
2n and 4n gridpoints, respectively. Without restriction of generality, we can assume that 
the gridspacing on a grid of n gridpoints is given hyh = 1/n. Assume furthermore that 
the discretized field ^„ is second order accurate, 

W„ = ^ + h^e2 + h^es + 0{h'^), (D.0.1) 

where the continuum error functions 62 and 63 do not depend on the gridspacing h. 
Then, 

W2„ = ^ + ih/2fe2 + ih/2)^e3 + (D.0.2) 

\iJ^„ = W + {h/4fe2 + ih/4fe3 + (D.0.3) 
We may subtract these gridfunctions on the common gridpoints, 

W2n - = -3/4h^e2 - IjSh^e^ + Oihl^), (D.0.4) 

and 

- = -3/4(///2)2e2 - 7/8ih/2fe4 + (D.0.5) 
Comparing the last two expressions we find that on the coarsest grid 

- ^2„ = 1/4 (1 + Om [W2„ - "¥„] . (D.0.6) 

By taking a discrete spatial norm over the grid, we can compute a time-dependent 
convergence factor 

0(t) ■= - ^2n\\ ,^ Q 7. 
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If the numerical scheme converges, then in the limit we should find 

limQ(0 = l/4. (D.0.8) 

In general we have, for h ^ 0, 

Wp2„ - Wpn = - ^„), (D.0.9) 

where a denotes the exponent of the leading error term in the discretized field ^„ = 
OQi") and |3 the factor between grid resolutions. 

A 2-level convergence test is applicable to discretized quantities that are supposed to 
be zero on the continuum level, i.e. 

W„ = + h^ez + h^e3 + 0{h% (D.0.10) 

It then follows that 

W2„ = 1/4.(1 +Oih))Wn. (D.0.11) 
More generally we have, for h ^ 0, 

^^jin = 1/n- (D.0.12) 
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